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Abstract 

Background: The primary visual cortex of many mammals contains a continuous representation of visual space, 
with a roughly repetitive aperiodic map of orientation preferences superimposed. It was recently found that 
orientation preference maps (OPMs) obey statistical laws which are apparently invariant among species widely 
separated in eutherian evolution. Here, we examine whether one of the most prominent models for the 
optimization of cortical maps, the elastic net (EN) model, can reproduce this common design. The EN model 
generates representations which optimally trade of stimulus space coverage and map continuity. While this model 
has been used in numerous studies, no analytical results about the precise layout of the predicted OPMs have 
been obtained so far. 

Results: We present a mathematical approach to analytically calculate the cortical representations predicted by the 
EN model for the joint mapping of stimulus position and orientation. We find that in all the previously studied 
regimes, predicted OPM layouts are perfectly periodic. An unbiased search through the EN parameter space 
identifies a novel regime of aperiodic OPMs with pinwheel densities lower than found in experiments. In an 
extreme limit, aperiodic OPMs quantitatively resembling experimental observations emerge. Stabilization of these 
layouts results from strong nonlocal interactions rather than from a coverage-continuity-compromise. 

Conclusions: Our results demonstrate that optimization models for stimulus representations dominated by 
nonlocal suppressive interactions are in principle capable of correctly predicting the common OPM design. They 
question that visual cortical feature representations can be explained by a coverage-continuity-compromise. 



Introduction 

The pattern of orientation columns in the primary visual 
cortex (VI) of carnivores, primates, and their close rela- 
tives are among the most intensely studied structures in 
the cerebral cortex and a large body of experimental (e. 
g., [1-13]) and theoretical work (e.g., [14-39]) has been 
dedicated to uncovering its organization principles and 
the circuit level mechanisms that underlie its develop- 
ment and operation. Orientation preference maps 
(OPMs) exhibit a roughly repetitive arrangement of pre- 
ferred orientations in which adjacent columns preferring 
the same orientation are separated by a typical distance 
in the millimeter range [2-5,10]. This range seems to be 
set by cortical mechanisms both intrinsic to a particular 
area [40] but potentially also involving interactions 
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between different cortical regions [41]. The pattern of 
orientation columns is however not strictly periodic 
because the precise local arrangement of preferred 
orientation never exactly repeats. Instead, OPMs appear 
as organized by a spatially complex aperiodic array of 
pinwheel centers, around which columns activated by 
different stimulus orientations are radially arranged like 
the spokes of a wheel [2-5,10]. The arrangement of 
these pinwheel centers, although spatially irregular, is 
statistically distinct from a pattern of randomly posi- 
tioned points [38] as well as from patterns of phase sin- 
gularities in a random pattern of preferred orientations 
[32,36,38,42] with spatial correlations identical to experi- 
mental observations [38,42]. This suggests that the lay- 
out of orientation columns and pinwheels although 
spatially aperiodic follows a definite system of layout 
rules. Cortical columns can in principle exhibit almost 
perfectly repetitive order as exemplified by ocular domi- 
nance (OD) bands in the macaque monkey primary 
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visual cortex [43,44]. It is thus a fundamental question 
for understanding visual cortical architecture, whether 
there are layout principles that prohibit a spatially 
exactly periodic organization of orientation columns and 
instead enforce complex arrangements of these columns. 

Recent comparative data have raised the urgency of 
answering this question and of dissecting what is consti- 
tutive of such complex layout principles. Kaschube et al. 
[38] quantitatively compared pinwheel arrangements in 
a large dataset from three species widely separated in 
the evolution of eutherian mammals. These authors 
found that the spatial statistics of pinwheels are surpris- 
ingly invariant. In particular, the overall pinwheel den- 
sity and the variability of pinwheel densities in regions 
from the scale of a single hypercolumn to substantial 
fractions of the entire primary visual cortex were found 
to be virtually identical. Characterizing pinwheel layout 
on the scale of individual hypercolumns, they found the 
distributions of nearest-neighbor pinwheel distances to 
be almost indistinguishable. Further supporting common 
layout rules for orientation columns in carnivores and 
primates, the spatial configuration of the superficial 
patch system [45] and the responses to drifting grating 
stimuli were recently found to be very similar in cat and 
macaque monkey primary visual cortex [46]. 

From an evolutionary perspective, the occurrence of 
quantitatively similar layouts for OPMs in primate tree 
shrews and carnivorous species appears highly informa- 
tive. The evolutionary lineages of these taxa diverged 
more than 65 million years ago during the basal radia- 
tion of eutherian mammals [47-49]. According to the 
fossil record and cladistic reconstructions, their last 
common ancestors (called the boreo-eutherial ancestors) 
were small-brained, nocturnal, squirrel-like animals of 
reduced visual abilities with a telencephalon containing 
only a minor neocortical fraction [47,50]. For instance, 
endocast analysis of a representative stem eutherian 
from the late cretaceous indicates a total anterior-pos- 
terior extent of 4 mm for its entire neocortex [47,50]. 
Similarly, the tenrec (Echinops telfari), one of the closest 
living relatives of the boreoeutherian ancestor [51,52], 
has a neocortex of essentially the same size and a visual 
cortex that totals only 2 mm 2 [47]. Since the neocortex 
of early mammals was subdivided into several cortical 
areas [47] and orientation hypercolumns measure 
between 0.4 and 1.4 mm 2 [38], it is difficult to envision 
ancestral eutherians with a system of orientation col- 
umns. In fact, no extant mammal with a visual cortex of 
such size is known to possess orientation columns [53]. 
It is therefore conceivable that systems of orientation 
columns independently evolved in laurasiatheria (such 
as carnivores) and in euarchonta (such as tree shrews 
and primates). Because galagos, tree-shrews, and ferrets 
strongly differ in habitat and ecologically relevant visual 



behaviors, it is not obvious that the quantitative similar- 
ity of pinwheel layout rules in their lineages evolved dri- 
ven by specific functional selection pressures (see [54] 
for an extended discussion). Kaschube et al. instead 
demonstrated that an independent emergence of identi- 
cal layout rules for pinwheels and orientation columns 
can be explained by mathematically universal properties 
of a wide class of models for neural circuit self- 
organization. 

According to the self-organization scenario, the com- 
mon design would result from developmental con- 
straints robustly imposed by adopting a particular kind 
of self-organization mechanism for constructing visual 
cortical circuitry. Even if this scenario is correct, one 
question still remains: What drove the different lineages 
to adopt a similar self-organization mechanism? As 
pointed out above, it is not easy to conceive that this 
adoption was favored by the specific demands of their 
particular visual habitats. It is, however, conceivable that 
general requirements for a versatile and representation- 
ally powerful cortical circuit architecture are realized by 
the common design. If this was true, the evolutionary 
benefit of meeting these requirements might have driven 
the adoption of large-scale self-organization and the 
emergence of the common design over evolutionary 
times. 

The most prominent candidate for such a general 
requirement is the hypothesis of a coverage-continuity- 
compromise (e.g., [19,21,55,56]). It states that the 
columnar organization is shaped to achieve an optimal 
tradeoff between the coverage of the space of visual sti- 
mulus features and the continuity of their cortical repre- 
sentation. On the one hand, each combination of 
stimulus features should be well represented in a corti- 
cal map to avoid 'blindness' to stimuli with particular 
feature combinations. On the other hand, the wiring 
cost to establish connections within the map of orienta- 
tion preference should be kept low. This can be 
achieved if neurons that are physically close in the cor- 
tex tend to have similar stimulus preferences. These two 
design goals generally compete with each other. The 
better a cortical representation covers the stimulus 
space, the more discontinuous it has to be. The tradeoff 
between the two aspects can be modeled in what is 
called a dimension reduction framework in which corti- 
cal maps are viewed as two-dimensional sheets which 
fold and twist in a higher-dimensional stimulus space 
(see Figure 1) to cover it as uniformly as possible while 
minimizing some measure of continuity [21,57,58]. 

From prior work, the coverage-continuity-compromise 
appears to be a promising candidate for a principle to 
explain visual cortical functional architecture. First, 
many studies have reported good qualitative agreement 
between the layout of numerically obtained dimension 
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Figure 1 The dimension reduction framework. In the dimension reduction framework, the visual cortex is modeled as a two-dimensional 
sheet that twists in a higher-dimensional stimulus (or feature) space to cover it as uniformly as possible while minimizing some measure of 
continuity (left). In this way, it represents a mapping from the cortical surface to the manifold of visual stimulus features such as orientation and 
retinotopy (right). 

L J 



reducing maps and experimental observations 
[19,21,42,55-57,59-71]. Second, geometric relationships 
between the representations of different visual features 
such as orientation, spatial frequency, and OD have 
been reproduced by dimension reduction models 
[25,56,63-65,67,68,72]. 

Mathematically, the dimension reduction hypothesis 
implies that the layouts of cortical maps can be under- 
stood as optima or near optima (global or local minima) 
of a free energy functional which penalizes both 'stimu- 
lus scotomas' and map discontinuity. Unfortunately, 
there is currently no dimension reduction model for 
which the precise layouts of optimal or nearly optimal 
solutions have analytically been established. To justify 
the conclusion that the tradeoff between coverage and 
continuity favors the common rules of OPM design 
found in experiment, knowledge of optimal dimension- 
reducing mappings however appears essential. 

Precise knowledge of the spatial organization of opti- 
mal and nearly optimal mappings is also critical for dis- 
tinguishing between optimization theories and frozen 
noise scenarios of visual cortical development. In a fro- 
zen noise scenario, essentially random factors such as 
haphazard wiring [73], the impact of spontaneous activ- 
ity patterns [74], or an idiosyncratic set of visual experi- 
ences [75] determine the emerging pattern of preferred 
orientations. This pattern is then assumed to be 'frozen' 
by an unknown mechanism, capable of preventing 
further modification of preferred orientations by 
ongoing synaptic turnover and activity-dependent plasti- 
city. Conceptually, a frozen noise scenario is diametri- 
cally opposed to any kind of optimization theory. Even 
if the reorganization of the pattern prior to freezing was 



to follow a gradient descent with respect to some cost 
function, the early stopping implies that the layout is 
neither a local nor a global minimum of this functional. 
Importantly, the layout of transient states is known to 
exhibit universal properties that can be completely inde- 
pendent of model details [25,32]. As a consequence, an 
infinite set of distinct optimization principles will gener- 
ate the same spatial structure of transient states. This 
implies in turn that the frozen transient layout is not 
specifically shaped by any particular optimization princi- 
ple. Map layouts will thus in principle only be informa- 
tive about design or optimization principles of cortical 
processing architectures if maps are not just frozen 
transients. 

In practice, however, the predictions of frozen noise 
and optimization theories might be hard to distinguish. 
Ambiguity between these mutually exclusive theories 
would result in particular, if the energy landscape of the 
optimization principle is so 'rugged' that there is essen- 
tially a local energy minimum next to any relevant ran- 
dom arrangement. Dimension reduction models are 
conceptually related to combinatorial optimization pro- 
blems like the traveling salesman problem (TSP) and 
many such problems are believed to exhibit rugged 
energy landscapes [76-78]. It is therefore essential to 
clarify whether paradigmatic dimension reduction mod- 
els are characterized by a rugged or a smooth energy 
landscape. If their energy landscapes were smooth with 
a small number of well-separated local minima, their 
predictions would be easy to distinguish from those of a 
frozen noise scenario. 

In this study, we examine the classical example of a 
dimension reduction model, the elastic network (EN) 
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model. Since the seminal work of Durbin and Mitchison 
[21], the EN model has widely been used to study visual 
cortical representations [25,42,62-65,69-72,79]. The EN 
model possesses an explicit energy functional which 
trades off a matching constraint which matches cortical 
cells to particular stimulus features via Hebbian learn- 
ing, with a continuity constraint that minimizes Eucli- 
dean differences in feature space between neighboring 
points in the cortex [63]. In two ways, the EN model's 
explicit variational structure is very appealing. First, cov- 
erage and continuity appear as separate terms in the 
free energy which facilitates the dissection of their rela- 
tive influences. Second, the free energy allows for the 
formulation of a gradient descent dynamics. The emer- 
gence of cortical selectivity patterns and their conver- 
gence toward a minimal energy state in this dynamics 
might serve as a model for an optimization process tak- 
ing place in postnatal development. 

Following Durbin and Mitchison, we consider the EN 
model for the joint mapping of two visual features: (i) 
position in visual space, represented in a retinotopic 
map (RM) and (ii) line orientation, represented in an 
OPM. To compute optimal dimension-reducing map- 
pings, we first develop an analytical framework for 
deriving closed-form expressions for fixed points, local 
minima, and optima of arbitrary optimization models 
for the spatial layout of OPMs and RMs in which pre- 
dicted maps emerge by a supercritical bifurcation as 
well as for analyzing their stability properties. By apply- 
ing this framework to different instantiations of the EN 
model, we systematically disentangle the effects of indi- 
vidual model features on the repertoire of optimal solu- 
tions. We start with the simplest possible model version, 
a fixed uniform retinotopy and an orientation stimulus 
ensemble with only a single orientation energy and then 
relax the uniform retinotopy assumption incorporating 
retinotopic distortions. An analysis for a second widely 
used orientation stimulus ensemble including also unor- 
iented stimuli is given in Appendix 1. Surprisingly, in all 
cases, our analysis yields pinwheel-free orientation 
stripes (OSs) or stereotypical square arrays of pinwheels 
as local minima or optimal orientation maps of the EN 
model. Numerical simulations of the EN confirm these 
findings. They indicate that more complex spatially 
aperiodic solutions are not dominant and that the 
energy landscape of the EN model is rather smooth. 
Our results demonstrate that while aperiodic stationary 
states exist, they are generally unstable in the considered 
model versions. 

To test whether the EN model is in principle capable 
of generating complex spatially aperiodic optimal orien- 
tation maps, we then perform a comprehensive unbiased 
search of the EN optima for arbitrary orientation stimu- 
lus distributions. We identify two key parameters 



determining pattern selection: (i) the intracortical inter- 
action range and (ii) the fourth moment of the orienta- 
tion stimulus distribution function. We derive complete 
phase diagrams summarizing pattern selection in the EN 
model for fixed as well as variable retinotopy. Small 
interaction ranges together with low to intermediate 
fourth moment values lead to pinwheel-free OSs, rhom- 
bic, or hexagonal crystalline orientation map layouts as 
optimal states. In the regime of large interaction ranges 
together with higher fourth moment, we find irregular 
aperiodic OPM layouts with low pinwheel densities as 
optima. Only in an extreme and previously unconsid- 
ered parameter regime of very large interaction ranges 
and stimulus ensemble distributions with an infinite 
fourth moment, optimal OPM layouts in the EN model 
resemble experimentally observed aperiodic pinwheel- 
rich layouts and quantitatively reproduce the recently 
described species-invariant pinwheel statistics. Unex- 
pectedly, we find that the stabilization of such layouts is 
not achieved by an optimal tradeoff between coverage 
and continuity of a localized population encoding by the 
maps but results from effectively suppressive long-range 
intracortical interactions in a spatially distributed repre- 
sentation of localized stimuli. 

We conclude our reexamination of the EN model with 
a comparison between different numerical schemes for 
the determination of optimal or nearly optimal map- 
pings. For large numbers of stimuli, numerically deter- 
mined solutions match our analytical predictions, 
irrespectively of the computational method used. 

Results and discussion 

Model definition and model symmetries 

We analyze the EN model for the joint optimization of 
position and orientation selectivity as originally intro- 
duced by Durbin and Mitchison [21]. In this model, the 
RM is represented by a mapping R(x) = C#i(x), R2M) 
which describes the receptive field center position of a 
neuron at cortical position x. Any RM can be decom- 
posed into an affine transformation x ^ X from cortical 
to visual field coordinates, on which a vector-field of 
retinotopic distortions r(x) is superimposed, i.e.: 

R(x) = X + r(x) 

with appropriately chosen units for x and R. 

The OPM is represented by a second complex-valued 
scalar field z(x). The pattern of orientation preferences 8 
(x) is encoded by the phase of z(x) via 

#( x ) = ^arg(z(x)). 

The absolute value \z(x)\ is a measure of the average 
cortical selectivity at position x. Solving the EN model 
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requires to find pairs of maps {r(x), z(x)} that represent 
an optimal compromise between stimulus coverage and 
map continuity. This is achieved by minimizing a free 
energy functional 

J=- = cr 2 C + K (1) 

in which the functional C measures the coverage of a 
stimulus space and the functional 71 the continuity of 
its cortical representation. The stimulus space is defined 
by an ensemble {S} of idealized point-like stimuli, each 
described by two features: s z = \s z \e 2ld and s r = (s x ,s y ) 
which specify the orientation 0 of the stimulus and its 
position in visual space s r (Figure 2b). C and 7Z are 
given by 

C[Z,T] = ~(ln f ^^-(| 5z -z( y )| 2 + |s r -X-r(y)|^)/2^| 
1 f 2 

j=i 

with V = (d XJ d y ) T , and 77 e [0, 1]. The ratios o 2 lr\ and 
a 2 lr\ r control the relative strength of the coverage term 
versus the continuity term for OPM and RM, respec- 
tively. (...) s denotes the average over the ensemble of 
stimuli. 

Minima of the energy functional T are stable fixed 
points of the gradient descent dynamics 

9 £ z(x) = - 2 S^=F*[z,r](x) 

5Z(X) (2) 

called the EN dynamics in the following. These 
dynamics read 

d t z(x) = ([s z - z(x)] e{x, S, z, r)) s + rj Az(x) (3) 

d t r(x) = ([s r - X - r(x)] e(x, S, z, r)) s + r) r Ar(x) , (4) 

where e(x, S, z, r) is the activity-pattern, evoked by a 
stimulus S = (s r , s z ) in a model cortex with retinotopic 
distortions r(x) and OPM z(x). It is given by 

^-(|s r -X-r(x)| 2 )/2a 2 ^-(|5 z -z(x)| 2 )/2a 2 
^ X "*'^ = / ^2 y ^-(|s r -X-r(y)| 2 )/2a^-(| 5z -z(y)| 2 )/2a2 ' 

Figure 2 illustrates the general features of the EN 
dynamics using the example of a single stimulus. Figure 
2a shows a model orientation map with a superimposed 
uniform representation of visual space. A single point- 
like, oriented stimulus S = (s r , s z ) with position s r = (s x , 
s y ) and orientation 6 - 1/2 arg(s z ) (Figure 2b) evokes a 



cortical activity pattern e(x, S, z, r) (Figure 2c). The 
activity-pattern in this example is of roughly Gaussian 
shape and is centered at the point, where the location s r 
of the stimulus is represented in cortical space. How- 
ever, depending on the model parameters and the sti- 
mulus, the cortical activity pattern may assume as well a 
more complex form (see also 'Discussion' section). 
According to Equations (3, 4), each stimulus and the 
evoked activity pattern induce a modification of the 
orientation map and RM, shown in Figure 2d. Orienta- 
tion preference in the activated regions is shifted toward 
the orientation of the stimulus. The representation of 
visual space in the activated regions is locally contracted 
toward the position of the stimulus. Modifications of 
cortical selectivities occur due to randomly chosen sti- 
muli and are set proportional to a very small learning 
rate. Substantial changes of cortical representations 
occur slowly through the cumulative effect of a large 
number of activity patterns and stimuli. These effective 
changes are described by the two deterministic equa- 
tions for the rearrangement of cortical selectivities equa- 
tions (3, 4) which are obtained by stimulus-averaging 
the modifications due to single activity patterns in the 
discrete stimulus model [25]. One thus expects that the 
optimal selectivity patterns and also the way in which 
cortical selectivities change over time are determined by 
the statistical properties of the stimulus ensemble. In 
the following, we assume that the stimulus ensemble 
satisfies three properties: (i) The stimulus locations s r 
are uniformly distributed across visual space, (ii) For the 
distribution of stimulus orientations, \s z \ and 6 are inde- 
pendent, (iii) Orientations 6 are distributed uniformly in 
[0, n\ 

These conditions are fulfilled by stimulus ensembles 
used in virtually all prior studies of dimension reduction 
models for visual cortical architecture (e.g., 
[19,21,25,64,65,71,72,80,81]). They imply several symme- 
tries of the model dynamics equations (3, 4). Due to the 
first property, the EN dynamics are equivariant under 
translations 

f y z(x) = z(x + y) 
f y r(x) = r(x + y), 

rotations 

Rpz(x) = e 2i ^z(Q^x) 
^r(x) = £V(^_/*x) 
with 2x2 rotation matrix 

/ cos£ -sin£\ 
fi cosfi )' 
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Figure 2 The EN model, (a) Example OPM (color code) together with a uniform map of visual space (RM) (grid lines), (b) Position s r = (s x , s y ) 
and orientation 0 of a 'pointlike' stimulus, (c) Cortical activity, evoked by the stimulus in b for the model maps in a. Dark regions are activated. 
Note, that in contrast to SOFM models, the activity pattern does not exhibit a stereotypical Gaussian shape, (d) Modification of orientation 
preference and retinotopy, caused by the stimulus in b. Orientation preferences prior to stimulus presentation are indicated with grey bars, after 
stimulus presentation with black bars. Most strongly modified preferences correspond to thick black bars. Modifications of orientation 
preferences and retinotopy are displayed on an exaggerated scale for illustration purposes. 



and reflections 



Pz(x) 
Pr(x) 



z(vl>x) 
^r(^x), 



where *F = diag(-l, 1) is the 2x2 reflection matrix. 
Equivariance means that 



f y F z [z,r] = F z [f y z,f y r] 



RpF z [z,v] = F[R pZt RpT] 



PF z [z,r] = F z [Pz,Pr], 



(5) 
(6) 
(7) 
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with mutatis mutandis the same equations fulfilled by 
the vector-field F r [z, r]. 

As a consequence, patterns that can be converted into 
one another by translation, rotation, or reflection of the 
cortical layers represent equivalent solutions of the 
model equations (3, 4), by construction. Due to the sec- 
ond assumption, the dynamics is also equivariant with 
respect to shifts in orientation S«pz(x) = e l<p z(x), i.e., 

e i4, F z [z,x] = F z [^z,r] (8) 

F r [z,r] = F r [^z,r]. (9) 

Thus, two patterns are also equivalent solutions of the 
model, if their layout of orientation domains and retino- 
topic distortions is identical, but the preferred orienta- 
tions differ everywhere by the same constant angle. 

Without loss of generality, we normalize the ensem- 
bles of orientation stimuli such that (|s z | 2 )s = {\s z \ 2 } = 2 
throughout this article. This normalization can always 
be restored by a rescaling of z(x) (see [25,69]). 

Our formulation of the dimension-reduction problem 
in the EN model utilizes a continuum description, both 
for cortical space and the set of visual stimuli. This facil- 
itates mathematical treatment and appears appropriate, 
given the high number of cortical neurons under one 
square millimeter of cortical surface (e.g., roughly 70000 
in cat VI [82]). Even an hypothesized neuronal mono- 
layer would consist of more than 20 x 20 neurons per 
hypercolumn area A 2 , constituting a quite dense sam- 
pling of the spatial periodicity. Treating the feature 
space as a continuum implements the concept that the 
cortical representation has to cover as good as possible 
the infinite multiplicity of conceivable stimulus feature 
combinations. 

The orientation unselective fixed point 

Two stationary solutions of the model can be estab- 
lished from symmetry. The simplest of these is the 
orientation unselective state with z(x) = 0 and uniform 
mapping of visual space r(x) = 0. First, by the shift sym- 
metry (Equation (8)), we find that z(x) = 0 is a fixed 
point of Equation (3). Second, by reflectional and rota- 
tional symmetry (Equations (5, 7)), we see that the 
right-hand side of Equation (4) has to vanish and hence 
the orientation unselective state with uniform mapping 
of visual space is a fixed point of Equations (3, 4). 

This homogeneous unselective state thus minimizes 
the EN energy functional if it is a stable solution of Equa- 
tions (3, 4). The stability can be determined by consider- 
ing the linearized dynamics of small deviations {r(x), z 
(x)} around this state. These linearized dynamics read 



9 t r(x) ~ Lr[r] = — ^ f d 2 ye ^ A r(y) + Vr A r(x) ( 10 ) 

167I(T 4 J 

3 t z(x) ~ L z [z] = (^2 -A *00 + r,Az{x) - —^4 / d 2 ye- ij ^-z{y), (11) 

where (A)# = {x t - y^Xj - yj) -2o% with <fy being Kro- 
necker's delta. We first note that the linearized 
dynamics of retinotopic distortions and orientation pre- 
ference decouple. Thus, up to linear order and near the 
homogeneous fixed point, both cortical representation 
evolve independently and the stability properties of the 
unselective state can be obtained by a separate examina- 
tion of the stability properties of both cortical 
representations. 

The eigenfunctions of the linearized retinotopy 
dynamics L r [r] can be calculated by Fourier-transform- 
ing Equation (10): 

2 

d t n(k) = -J2 (^e-^hkj + rjrk^ij) ?j(k), 
i=i 

where k = |k| and i = 1, 2. A diagonalization of this 
matrix equation yields the eigenvalues 

X[ = -k 2 {rj r + e- alkl o 2 ), X r T = -rjrk 2 
with corresponding eigenfunctions (in real space) 

r L (x) = V'^ + cc. 
r T (x) = k 0W2 A x + ex., 

where k«p = |k|(cos (p, sin (p) T . These eigenfunctions 
are longitudinal (L) or transversal (T) wave patterns. In 
the longitudinal wave, the retinotopic distortion vector r 
(x) lies parallel to k which leads to a 'compression wave' 
(Figure 3b, left). In the transversal wave pattern (Figure 
3b, right), the retinotopic distortion vector is orthogonal 
to k. We note that the linearized Kohonen model [61] 
was previously found to exhibit the same set of eigen- 
functions [80]. Because both spectra of eigenvalues U T , 
A^are smaller than zero for every a >0, r\ r >0, and k >0 
(Figure 3a), the uniform retinotopy r(x) = 0 is a stable 
fixed point of Equation (4) irrespective of parameter 
choice. 

The eigenfunctions of the linearized OPM dynamics 
L z [z] are Fourier modes ~ e lk * by translational symme- 
try. By rotational symmetry, their eigenvalues only 
depend on the wave number k and are given by 

X z (k) = -1 + (l - <r feV ) - r]k 2 
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Figure 3 The linearization of the EN model dynamics around the unselective fixed point, (a) Eigenvalue spectra of the linearized 
retinotopy dynamics for longitudinal mode {k r L (k) , blue trace) and transversal mode {X r T (k), red trace), (b) Longitudinal mode 
~ k00 ik ^ x + ex. (left) and transversal mode ~ k0 +7r /2^ lk</>x + CC (right), (c) Spectrum of eigenvalues of the linearized OPM dynamics (red 
trace) for a < o*(ry). Orange region marks the unstable annulus of Fourier modes (critical circle), (d) Stability regions of the nonselective state in 
the EN model. The stability border is given by 0^(77) (Equation (12)). 



(see [25]). This spectrum of eigenvalues is depicted in 
Figure 3c. For 77 >0, X Z (I<) has a single maximum at 

k c = ^ln(l/r,). For 

a > cr*[rj) = + rjlnrj — rj , (12) 

this maximal eigenvalue r - A, z (k c ) is negative. Hence, 
the unselective state with uniform retinotopy is a stable 
fixed point of Equations (3, 4) and the only known solu- 
tion of the EN model in this parameter range. For a < 
c7*(77), the maximal eigenvalue r is positive, and the non- 
selective state is unstable with respect to a band of 
Fourier modes ~ e lkx with wave numbers around |k| ~ 
k c (see Figure 3c). This annulus of unstable Fourier 
modes is called the critical circle. The finite wavelength 
instability [83-85] (or Turing instability [86]) leads to 
the emergence of a pattern of orientation preference 
with characteristic spacing A = 2n/k c from the nonselec- 
tive state on a characteristic timescale r = II r. 



One should note that as in other models for the self- 
organization of orientation columns, e.g., [15,57], the 
characteristic spatial scale A arises from effective intra- 
cortical interactions of 'Mexican-hat' structure (short- 
range facilitation, longer-ranged suppression). The 
short-range facilitation in the linearized EN dynamics is 
represented by the first two terms on the right-hand 
side of Equation (11). Since a <1 in the pattern forming 
regime, the prefactor in front of the first term is posi- 
tive. Due to the second, Laplacian term, it is favored 
that neighboring units share selectivity properties, a pro- 
cess mediated by short-range facilitation. Longer-ranged 
suppression is represented by the convolution term in 
Equation (11). 

Mathematically, this term directly results from the 
soft-competition in the 'activity-dependent' coverage 
term of Equation (1). The local facilitation is jointly 
mediated by coverage (first term) and continuity (second 
term) contributions. 
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Figure 3d summarizes the result of the linear stability 
analysis of the nonselective state. For a > (7*(77), the 
orientation unselective state with uniform retinotopy is 
a minimum of the EN-free energy and also the global 
minimum. For a < <7*(/7), this state represents a maxi- 
mum of the energy functional and the minima must 
thus exhibit a space-dependent pattern of orientation 
selectivities. 

Orientation stripes 

Within the potentially infinite set of orientation selective 
fixed points of the model, one class of solutions can be 
established from symmetry: {r(x) = 0, z(x) = A 0 e lklc }. In 
these pinwheel-free states, orientation preference is con- 
stant along one axis in cortex (perpendicular to the vec- 
tor k), and each orientation is represented in equal 
proportion (see Figure 4a). Retinotopy is perfectly uni- 
form. Although this state may appear too simple to be 
biologically relevant, we will see that it plays a funda- 
mental role in the state space of the EN model. It is 
therefore useful to establish its existence and basic char- 
acteristics. The existence of OS solutions follows directly 
from the model's symmetries (Equations (5) to (9)). 
Computing 

0]] = F z [T y [e ikx ], T y [0]] = F z [^ ikx , 0] = e iky F z [e ikx / 0] 

demonstrates that i^[e lkx , 0] is proportional to e' kx . 
This establishes that the subspace of functions ~ e' kx is 
invariant under the dynamics given by Equation (3). For 
A 0 = 0, we recover the trivial fixed point of the EN 
dynamics by construction, as shown above. This means 
that within this subspace A 0 = 0 is either a minimum or 
a maximum of the EN energy functional (Equation (1)). 
Furthermore, for A 0 — > ©o the EN energy tends to infi- 
nity. If the trivial fixed point is unstable, it corresponds 
to a maximum of the EN energy functional. Therefore, 
there must exist at least one minimum with A 0 * 0 in 
the subspace of functions ~ e zkx which then corresponds 
to a stationary state of the EN dynamics. 

Regarding the dynamics of retinotopic deviations, the 
model's symmetries equations can be invoked to show 
that for the state {0, A 0 e /kx }, the right-hand side of 
Equation (4) has to be constant in space: 

T y [F r > lkx / 0]] = F r [T y [g lkx ],T y [0]] = F r [^ lkx ,0] = F^^O] 

If this constant was nonzero the RM would drift with 
constant velocity. This, however, is impossible in a var- 
iational dynamics such that this constant must vanish. 
The OS solution (Figure 4a) is to the best of our knowl- 
edge the only exact nontrivial stationary solution of 
Equations (3, 4) that can be established without any 
approximations. 



Doubly periodic and quasi-periodic solutions 

In the EN model as considered in this study, the maps 
of visual space and orientation preference are jointly 
optimized to trade off coverage and continuity leading 
to mutual interactions between the two cortical repre- 
sentations. These mutual interactions vanish in the rigid 
retinotopy limit T] r — > ©o and the perfectly uniform reti- 
notopy becomes an optimal solution for arbitrary orien- 
tation column layout z(x). As it is not clear how 
essential the mutual interactions with position specificity 
are in shaping the optimal orientation column layout, 
we continue our investigation of solution classes by con- 
sidering global minima of optimization models with 
fixed uniform retinotopy. The mutual interactions will 
be taken into account in a subsequent step. 

In the rigid retinotopy limit, minima of the energy 
functional are stable stationary states of the dynamics of 
the OPM (Equation (3)) with r(x) = 0. To compute 
orientation selective stationary solutions of this OPM 
dynamics, we employ that in the vicinity of a supercriti- 
cal bifurcation where the nonselective fixed point z(x) = 
0 becomes unstable, the entire set of nontrivial fixed 
points is determined by the third-order terms of the 
Volterra series representation of the operator F*[z, 0] 
[35,84,85,87]. The symmetries given by Equations (5) to 
(9) restrict the general form of such a third-order 
approximation for any model of OPM optimization to 



d t z(x) ^ L z [z] + N z 3 [z, z, z], 



(13) 



where the cubic operator N 3 is written in trilinear 
form, i.e., 



a i z i' P kZk ' Y YlZl 



= ^ajPkyiNllzpZktZi] 



In particular, all even terms in the Volterra Series 
representation of ^[z, 0] vanish due to the Shift-Sym- 
metry (Equations (8, 9)). Explicit analytic computation 
of the cubic nonlinearities for the EN model is cumber- 
some but not difficult (see 'Methods' section) and yields 
a sum 



N3 [z, z, z] = djN ] 3 [z, z, z] . 



(14) 



The individual nonlinear operators M ] 3 are with one 

exception nonlocal convolution-type operators and are 
given in the 'Methods' section (Equation (38)), together 
with a detailed description of their derivation. 

Only the coefficients <Zj depend on the properties of 
the ensemble of oriented stimuli. 
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Figure 4 Exact and approximate orientation selective fixed points of OPM optimization models, (a) Pinwheel-free OS pattern. Diagram 
shows the position of the wave vector in Fourier space, (b) rPWC with four nonzero wave vectors, (c) Essentially complex planforms (ECPs). The 
index n indicates the number of nonzero wave vectors. The index / enumerates nonequivalent configurations of wave vectors with the same n, 
starting with /' = 0 for the most anisotropic planform. For n = 3, 5, and 15, there are 2, 4, and 612 different ECPs, respectively. OPM layouts 
become more irregular with increasing n. 
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To calculate the fixed points of Equation (13), we use 
a perturbative method called weakly nonlinear analysis 
that enables us to analytically examine the structure and 
stability of inhomogeneous stationary solutions in the 
vicinity of a finite-wavelength instability. Here, we exam- 
ine the stability of so-called planforms [83-85]. Plan- 
forms are patterns that are composed of a finite number 
of Fourier components, such as 



ik,x 



z(x) = £>(ty 



for a pattern of orientation columns. With the above 
planform ansatz, we neglect any spatial dependency of 
the amplitudes A ; (t) for example due to long-wave 
deformations for the sake of simplicity and analytical 
tractability. When the dynamics is close to a finiite 
wavelength instability, the essential Fourier components 
of the emerging pattern are located on the critical circle 
|k ; | = k c . The dynamic equations for the amplitudes of 
these Fourier components are called amplitude equa- 
tions. For a discrete number of N Fourier components 
of z(x) whose wave vectors lie equally spaced on the cri- 
tical circle, the most general system of amplitude equa- 
tions compatible with the model's symmetries 
(Equations (5) to (9)) has the form [35,87] 



M = rA x - A{ J2gij\Aj\ 2 - EMA-' 



(15) 



with r >0. Here, g t j and^y are the real-valued coupling 
coefficients between the amplitudes A t and Aj. They 
depend on the differences between indices \i - j\ and are 
entirely determined by the nonlinearity N3[z,z,z] in 
Equation (13). If the wave vectors = (cos a b sin a-)k c 
are parameterized by the angles a b then the coefficients 
gij andfij are functions only of the angle a - \a t - aj\ 
between the wave vectors k/ and k ; . One can thus obtain 
the coupling coefficients from two continuous functions 
g(a) and f(a) that can be obtained from the nonlinearity 
N3[z,z,z] (see 'Methods' section for details). In the fol- 
lowing, these functions are called angle-dependent inter- 
action functions. The amplitude equations are 
variational if and only if g t j and/y are real-valued. In this 
case they can be derived through 

A dAj 
from an energy 



If the coefficients g t j and are derived from Equation 
(1), the energy U A for a given planform solution corre- 
sponds to the energy density of the EN energy func- 
tional considering only terms up to fourth-order in z(x). 

The amplitude equations (15) enable to calculate an 
infinite set of orientation selective fixed points. For the 
above OS solution with one nonzero wave vector z(x) = 
A 0 e' kx , the amplitude equations predict the so far unde- 
termined amplitude 



|A 0 | 



and its energy 



Uqs 



2gu 



(17) 



(18) 



Since g u >0, this shows that OS stationary solutions 
only exist for r >0, i.e., in the symmetry breaking 
regime. As for all following fixed-points, U os specifies 
the energy difference to the homogeneous unselective 
state z(x) = 0. 

A second class of stationary solutions can be found 
with the ansatz 



z(x) 



: Ax^' klX + A 2 ^' k2X + A 3 £T iklX + A 4 e~ ik2X 



with amplitudes Aj = \Aj\e l(Pj and ^(ki, k 2 ) = ci >0. By 
inserting this ansatz into Equation (15) and assuming 
uniform amplitude |Ai| = |A 2 | = |A 2 | = |A 4 | =A, 
we obtain 



A 2 



gOO + gOn + gOa + gOiz-a ~ 2foa 



(19) 



The phase relations of the four amplitudes are given 



by 



h + 93 

b 2 + 04 



00 

00 + X. 



These solutions describe a regular rhombic lattice of 
pinwheels and are therefore called rhombic pinwheel 
crystals (rPWCs) in the following. Three phases can be 
chosen arbitrarily according to the two above condi- 
tions, e.g., (p 0 , A 0 = (pi - (p 3 and A x = q> 2 - (p4< For an 
rPWC parameterized by these phases, A 0 shifts the 
absolute positions of the pinwheels in ^-direction, A x 
shifts the absolute positions of the pinwheels in j-direc- 
tion, and q> 0 shifts all the preferred orientations by a 
constant angle. The energy of an rPWC solution is 



U A = -rjr\ Al f + \ ^lAdW + I jhfijAiAi-AjAj-. (16) U ™ C = 



2r 



i,j=l 



i,j=l 



gOO + gOn + gOa + gOn -a 



Oa 



(20) 
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An example of such a solution is depicted in Figure 
4b. We note that rPWCs have been previously found in 
several other models for OPM development 
[27,31,37,39,88]. The pinwheel density p of an rPWC, i. 
e., the number of pinwheels in an area of size A 2 , is 
equal to p = 4 sin a [54]. The angle a which minimizes 
the energy U rPWC can be computed by maximizing the 
function 

= gOa + gOn-a ~ 2foa (21) 

in the denominator of Equation (20). 

The two solution classes discussed so far, namely OS 
and rPWCs, exhibit one prominent feature, absent in 
experimentally observed cortical OPMs, namely perfect 
spatial periodicity. Many cortical maps including OPMs 
do not resemble a crystal-like grid of repeating units. 
Rather the maps are characterized by roughly repetitive 
but aperiodic spatial arrangement of feature preferences 
(e.g., [5,10]). This does not imply that the precise layout 
of columns is arbitrary. It rather means that the rules of 
column design cannot be exhaustively characterized by 
mapping a representative' hypercolumn. 

Previous studies of abstract models of OPM develop- 
ment introduced the family of so-called essentially com- 
plex planforms (ECPs) as stationary solutions of 
Equation (15). This solution class encompasses a large 
variety of realistic quasi-periodic OPM layouts and is 
therefore a good candidate solution class for models of 
OPM layouts. In addition, Kaschube et al. [38] demon- 
strated that models in which these are optimal solutions 
can reproduce all essential features of the common 
OPM design in ferret, tree-shrew, and galago. An n-ECP 
solution can be written as 

n 

Z (x) = ^ A/^ X , 
i=i 

with n - N/2 wave vectors k y = k c (cos(n;j/n), sm(jij/n)) 
distributed equidistantly on the upper half of the critical 
circle, complex amplitudes Aj and binary variables lj = 
±1 determining whether the mode with wave vector k ; 
or -k ; is active (nonzero). Because these planforms can- 
not realize a real-valued function they are called essen- 
tially complex [35]. For an n-ECP, the third term on the 
right-hand side of Equation (15) vanishes and the ampli- 
tude equations for the active modes A t reduce to a sys- 
tem of Landau equations 

n 

A i = rA i -A l Y J gim 2 , 

i=i 

where gy is the n x n-coupling matrix for the active 
modes. Consequently, the stationary amplitudes obey 



n 

Wf-'Efe" 1 )*' (22) 
The energy of an n-ECP is given by 
Uecp = - ^E(3 -1 )if (23) 

a 

We note that this energy in general depends on the 
configuration of active modes, given by the //s, and 
therefore planforms with the same number of active 
modes may not be energetically degenerate. 

Families of n-ECP solutions are depicted in Figure 4c. 
The 1-ECP corresponds to the pinwheel-free OS pattern 
discussed above. For fixed n > 3, there are multiple 
planforms not related by symmetry operations which 
considerably differ in their spatial layouts. For n > 4, the 
patterns are spatially quasi-periodic, and are a generali- 
zation of the so-called Newell-Pomeau turbulent crystal 
[89,90]. For n > 10, their layouts resemble experimen- 
tally observed OPMs. Different n-ECPs however differ 
considerably in their pinwheel density. Planforms whose 
nonzero wave vectors are distributed isotropically on the 
critical circle typically have a high pinwheel density (see 
Figure 4c, n = 15 lower right). Anisotropic planforms 
generally contain considerably fewer pinwheels (see Fig- 
ure 4c, n = 15 lower left). All large n-ECPs, however, 
exhibit a complex quasi-periodic spatial layout and a 
nonzero density of pinwheels. 

In order to demonstrate that a certain planform is an 
optimal solution of an optimization model for OPM lay- 
outs in which patterns emerge via a supercritical bifur- 
cation, we not only have to show that it is a stationary 
solution of the amplitude equations but have to analyze 
its stability properties with respect to the gradient des- 
cent dynamics as well as its energy compared to all 
other candidate solutions. 

Many stability properties can be characterized by 
examining the amplitude equations (15). In principle, 
the stability range of an n-ECPs may be bounded by two 
different instability mechanisms: (i) an intrinsic instabil- 
ity by which stationary solutions with n active modes 
decay into ones with lower n. (ii) an extrinsic instability 
by which stationary solutions with a 'too low' number of 
modes are unstable to the growth of additional active 
modes. These instabilities can constrain the range of 
stable n to a small finite set around a typical n [35,87]. 
A mathematical evaluation of both criteria leads to pre- 
cise conditions for extrinsic and intrinsic stability of a 
planform (see 'Methods' section). In the following, a 
planform is said to be stable, if it is both extrinsically 
and intrinsically stable. A planform is said to be an opti- 
mum (or optimal solution) if it is stable and possesses 
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the minimal energy among all other stationary planform 
solutions. 

Taken together, this amplitude equation approach 
enables to analytically compute the fixed points and 
optima of arbitrary optimization models for visual corti- 
cal map layout in which the functional architecture is 
completely specified by the pattern of orientation col- 
umns z(x) and emerges via a supercritical bifurcation. 
Via a third-order expansion of the energy functional 
together with weakly nonlinear analysis, the otherwise 
analytically intractable partial integro-differential equa- 
tion for OPM layouts reduces to a much simpler system 
of ordinary differential equations, the amplitude equa- 
tions. Using these, several families of solutions, OSs, 
rPWCs, and essentially complex planforms, can be sys- 
tematically evaluated and comprehensively compared to 
identify sets of unstable, stable and optimal, i.e., lowest 
energy fixed points. As already mentioned, the above 
approach is suitable for arbitrary optimization models 
for visual cortical map layout in which the functional 
architecture is completely specified by the pattern of 
orientation columns z(x) which in the EN model is ful- 
filled in the rigid retinotopy limit. We now start by con- 
sidering the EN optimal solutions in this limit and 
subsequently generalize this approach to models in 
which the visual cortical architecture is jointly specified 
by maps of orientation and position preference that are 
matched to one another. 

Representing an ensemble of 'bar'-stimuli 

We start our investigation of optimal dimension-redu- 
cing mappings in the EN model using the simplest and 
most frequently used orientation stimulus ensemble, the 
distribution with s^-values uniformly arranged on a ring 
with radius r Sz = \fl [57, 64-66,91]. We call this stimulus 
ensemble the circular stimulus ensemble in the follow- 
ing. According to the linear stability analysis of the non- 
selective fixed point, at the point of instability, we 
choose a = (7* (77) such that the linearization given in 
Equation (11) is completely characterized by the conti- 
nuity parameter 77. Equivalent to specifying 77 is to fix 
the ratio of activation range a and column spacing A 

a/A = J-,/log(l/i7) (24) 

as a more intuitive parameter. This ratio measures the 
effective interaction-range relative to the expected spa- 
cing of the orientation preference pattern. In abstract 
optimization models for OPM development a similar 
quantity has been demonstrated to be a crucial determi- 
nant of pattern selection [35,87]. We note, however, that 
due to the logarithmic dependence of o7A on 77, a slight 
variation of the effective interaction range may 



correspond to a variation of the continuity parameter 77 
over several orders of magnitude. In order to investigate 
the stability of stationary planform solutions in the EN 
model with a circular orientation stimulus ensemble, we 
have to determine the angle-dependent interaction func- 
tions g{a) and f(a). For the coefficients aj in Equation 
(14) we obtain 

a _ _J L 4. _L_ n - -J L_ n - 1 , 1 

U\ 4o . 6 4 + 2a2 U2 4 5 g 8 "3 167TCT 8 + 87TCT 6 

n = L_ . _J L_ n = 1 n _ _J i__ 

" 4 8na 8 + 47TO-6 Sna A " 5 16ttct 8 U<s 8ttct 6 g 16ttct 8 

d? = I2n 2 a w ~ I2n 2 a 8 ^ 8 = 24tt 2 ct 10 ^ = ~~ 64tt 3 ct 12 

dl ° = 12tt 2 ct 10 _ 12n 2 a 8 dn = 2An 2 o w ' 

The angle-dependent interaction functions of the EN 
model with a circular orientation stimulus ensemble are 
then given by 

g{a) - J ? (l-2^-^(-- 1 )(l-2^ 2 — )) 

+ (^(cos^) _ ^ + ^ W sinh4(1/2 ^ C osa) 

f{a) = ^ (l -e- 2k '° 2 (cosh(2^ (T 2 cosa) +2cosh(fe c 2 CT 2 cosQ ; )) + 2g- fe ' CT2 ) 

+ (e- 2k ^ 2 cosh(2fe 2 a 2 cos a) - l) + ^~ lkW sinh 4 (l/2fe c 2 a 2 cos a) . 

These functions are depicted in Figure 5 for two dif- 
ferent values of the interaction range o7A. We note that 
both functions are positive for all a/ A which is a suffi- 
cient condition for a supercritical bifurcation from the 
homogeneous nonselective state in the EN model. 

Finally, by minimizing the function s(a) in Equation 
(21), we find that the angle a which minimizes the 
energy of the rPWC fixed-point is a = tt/2. This corre- 
sponds to a square array of pinwheels (sPWC). Due to 
the orthogonal arrangement oblique and cardinal orien- 
tation columns and the maximized pinwheel density of 
p = 4, the square array of pinwheels has the maximal 
coverage among all rPWC solutions. 
Optimal solutions close to the pattern formation threshold 
We first tested for the stability of pinwheel-free OS 
solutions and the sPWCs, by analytical evaluation of the 
criteria for intrinsic and extrinsic stability (see 'Methods' 
section). We found both, OSs and sPWCs, to be intrinsi- 
cally and extrinsically stable for all o7A. Next, we tested 
for the stability of n-ECP solutions with 2 < n < 20. We 
found all w-ECP configurations with 2 < n < 20 to be 
intrinsically unstable for all cr/A. Hence, none of these 
planforms represent optimal solutions of the EN model 
with a circular stimulus ensemble, while both OSs and 
sPWC are always local minima of the energy functional. 

By evaluating the energy assigned to the sPWC (Equa- 
tion 20) and the OS pattern (Equation 18), we next 
identified two different regimes: (i) For short interaction 
range a I A ^ 0.122 the sPWC possesses minimal energy 
and is therefore the predicted global minimum, (ii) For 
a/A > 0.122 the OS pattern is optimal. 

Figure 6a shows the resulting simple phase diagram. 
sPWCs and OSs are separated by a phase border at a/A 
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Figure 6 Optimal solutions of the EN model with a circular orientation stimulus ensemble [57,64-66,91 ]and fixed representation of 
visual space, (a) At criticality, the phase space of this model is parameterized by either the continuity parameter r\ (blue labels) or the 
interaction range cr/A (red labels, see text), (b, c) OPMs (b) and their power spectra (c) in a simulation of Equation (3) with r(x) = 0 and r = 0.1, 
o/A = 0.1 (77 = 0.67) and circular stimulus ensemble (see also Additional file 1). (d) Analytically predicted optimum for cr/A < 0.122 (quadratic 
pinwheel crystal), (e) Pinwheel density time courses for four different simulations (parameters as in b; gray traces, individual realizations; black 
trace, simulation in b; red trace, mean value), (f) Mean squared amplitude of the stationary pattern, obtained in simulations (parameters as in b) 
for different values of the control parameter r (black circles) and analytically predicted value (solid green line), (g, h) OPMs (g) and their power 
spectra (h) in a simulation of Equation (3) with r(x) = 0 and o/A = 0.15 (77 = 0.41) (other parameters as in b, see also Additional file 2). (i) 
Analytically predicted optimum for cr/A > 0.122 (orientation stripes), (j) Pinwheel density time courses for four different simulations (parameters 
as in g; gray traces, individual realizations; black trace, simulation in g; red trace, mean value), (k) Mean squared amplitude of the stationary 
pattern, obtained in simulations (parameters as in g) for different values of the control parameter r (black circles) and analytically predicted value 
(solid green line). 
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« 0.122. We numerically confirmed these analytical pre- 
dictions by extensive simulations of Equation (3) with r 
(x) = 0 and the circular stimulus ensemble (see 'Meth- 
ods' section for details). Figure 6b,c shows snapshots of 
a representative simulation with short interaction range 
(r = 0.1, a/A = 0.1 (77 = 0.67)) (see also Additional file 
1). After the phase of initial pattern emergence (symme- 
try breaking), the OPM layout rapidly approaches a 
square array of pinwheels, the analytically predicted 
optimum (Figure 6d). Pinwheel density time courses 
(see 'Methods' section) display a rapid convergence to a 
value close to the predicted density of 4 (Figure 6e). Fig- 
ure 6f shows the stationary mean squared amplitudes of 
the pattern obtained for different values of the control 
parameter r (black circles). For small control para- 
meters, the pattern amplitude is perfectly predicted by 
Equation (19) (solid green line). Figure 6g,h shows snap- 
shots of a typical simulation with longer interaction 
range (r = 0.1, a/A = 0.15 (rj = 0.41)) (see also Addi- 
tional file 2). After the emergence of an OPM with 
numerous pinwheels, pinwheels undergo pairwise anni- 
hilation as previously described for various models of 
OPM development and optimization [25,27,35]. The OP 
pattern converges to a pinwheel-free stripe pattern, 
which is the analytically computed optimal solution in 
this parameter regime (Figure 6i). Pinwheel densities 
decay toward zero over the time course of the simula- 
tions (Figure 6j). Also in this parameter regime, the 
mean squared amplitude of the pattern is well-predicted 
by Equation (17) for small r (Figure 6k). 

In summary, the phase diagram of the EN model with 
a circular stimulus ensemble close to threshold is 
divided into two regions: (i) for a small interaction 
range (large continuity parameter) a square array of pin- 
wheels is the optimal dimension-reducing mapping and 
(ii) for a larger interaction range (small continuity para- 
meter) OSs are the optimal dimension-reducing map- 
ping. Both states are stable throughout the entire 
parameter range. All other planforms, in particular 
quasi-periodic n-ECPs are unstable. At first sight, this 
structure of the EN phase diagram may appear rather 
counterintuitive. A solution with many pinwheel-defects 
is energetically favored over a solution with no defects 
in a regime with large continuity parameter where dis- 
continuity should be strongly penalized in the EN 
energy term. However, a large continuity parameter at 
pattern formation threshold inevitably leads to a short 
interaction range a compared to the characteristic spa- 
cing A (see Equation (24)). In such a regime, the gain in 
coverage by representing many orientation stimuli in a 
small area spanning the typical interaction range, e.g., 
with a pinwheel, is very high. Our results show that the 
gain in coverage by a spatially regular positioning of 



pinwheels outweighs the accompanied loss in continuity 
above a certain value of the continuity parameter. 
EN dynamics far from pattern formation threshold 
Close to pattern formation threshold, we found only two 
stable solutions, namely OSs and sPWCs. Neither of the 
two exhibits the characteristic aperiodic and pinwheel- 
rich organization of experimentally observed OPMs. 
Furthermore, the pinwheel densities of both solutions (p 
= 0 for OSs and p = 4 for sPWCs) differ considerably 
from experimentally observed values [38] around 3.14. 
One way toward more realistic stable stationary states 
might be to increase the distance from pattern forma- 
tion threshold. In fact, further away from threshold, our 
perturbative calculations may fail to correctly predict 
optimal solutions of the model due to the increasing 
influence of higher order terms in the Volterra series 
expansion of the right-hand side in Equation (3). 

To asses this possibility, we simulated Equation (3) 
with r(x) = 0 and a circular stimulus ensemble for very 
large values of the control parameter r. Figure 7 dis- 
plays snapshots of such a simulation for r = 0.8 as well 
as their pinwheel density time courses for two different 
values of o7A. Pinwheel annihilation in the case of 
large o7A is less rapid than close to threshold (Figure 
7a,b). The OPM nevertheless converges toward a lay- 
out with rather low pinwheel density with pinwheel- 
free stripe-like domains of different directions joined 
by domains with essentially rhombic crystalline pin- 
wheel arrangement. The linear zones increase their 
size over the time course of the simulations, eventually 
leading to stripe-patterns for large simulation times. 
For smaller interaction ranges o7A, the OPM layout 
rapidly converges toward a crystal-like rhombic 
arrangement of pinwheels, however containing several 
dislocations (Figure 24a in Appendix 1) [84]. Disloca- 
tions are defects of roll or square patterns, where two 
rolls or squares merge into one, thus increasing the 
local wavelength of the pattern [83,85]. Nevertheless, 
for all simulations, the pinwheel density rapidly 
reaches a value close to 4 (Figure 7c) and the square 
arrangement of pinwheels is readily recognizable. Both 
features, the dislocations in the rhombic patterns and 
domain walls in the stripe patterns, have been fre- 
quently observed in pattern-forming systems far from 
threshold [84,85]. 

In summary, the behavior of the EN dynamics with 
circular stimulus ensemble far from pattern formation 
threshold agrees very well with our analytical predic- 
tions close to threshold. Again, orientation stripes and 
square pinwheel crystals are identified as the only sta- 
tionary solutions. Aperiodic and pinwheel-rich patterns 
which resemble experimentally OPM layouts were not 
observed. 
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Figure 7 Numerical analysis of the EN dynamics with circular orientation stimulus ensemble and fixed representation of visual space 
far from pattern formation threshold, (a) OPMs and their power spectra in a simulation of Equation (3) with r(x) = 0, r = 0.8, a/A = 0.3 (77 = 

0.028) and circular orientation stimulus ensemble. Pinwheel density time courses for four different simulations (parameters as in a; gray traces, 

individual realizations; black trace, simulation in a; red trace, mean value) (c, d) OPMs and their power spectra in a simulation of Equation (3) 

with r(x) = 0, r = 0.8, a/A = 0.12 (77 = 0.57) and circular orientation stimulus ensemble, (d) Pinwheel density time courses for four different 

simulations (parameters as in c; gray traces, individual realizations; black trace, simulation in c; red trace, mean value). 
^ J 



Taking retinotopic distortions into account 

So far, we have examined the optimal solutions of the 
EN model for the simplest and most widely used orien- 
tation stimulus ensemble. Somewhat unexpected from 
previous reports, the optimal states in this case do not 
exhibit the irregular structure of experimentally 
observed orientation maps. Our treatment however dif- 
fers from previous approaches in that the mapping of 
visual space so far was assumed to be undistorted and 
fixed, i.e., r(x) = 0. We recall that in their seminal publi- 
cation, Durbin and Mitchison [21] in particular demon- 
strated interesting correlations between the map of 
orientation preference and the map of visual space. 
These correlations suggest a strong coupling between 
the two that may completely alter the model's dynamics 
and optimal solutions. 

It is thus essential to clarify whether the behavior of 
the EN model observed above changes or persists if we 
relax the simplifying assumption of undistorted retino- 
topy and allow for retinotopic distortions. By analyzing 
the complete EN model dynamics (Equations (3, 4)), we 



study the EN model exactly as originally introduced by 
Durbin and Mitchison [21]. 

We again employ the fact that in the vicinity of a 
supercritical bifurcation where the nonorientation selec- 
tive state becomes unstable, the entire set of nontrivial 
fixed points of Equations (3, 4) is determined by the 
third-order terms of the Volterra series representation 
of the nonlinear operators F*[z, r] and F r [z, r]. The 
model symmetries equations (5) to (9) restrict the gen- 
eral form of the leading order terms for any model for 
the joint optimization of OPM and RM to 

d t z(x) = L z [z] + Q z [r, z] + N3 [z, z, z] + • • • (26) 

3 t r(x) = L r [r] + Q r [z,z] + --. . (27) 

Because the uniform retinotopy is linearly stable, reti- 
notopic distortions are exclusively induced by a coupling 
of the RM to the OPM via the quadratic vector-valued 
operator Q r [z,z]. These retinotopic distortions will in 
turn alter the dynamics of the OPM via the quadratic 
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complex-valued operator Q z [r, z\. Close to the point of 
pattern onset (r « 1), the timescale of OPM develop- 
ment, i = II r, becomes arbitrarily large and retinotopic 
deviations evolve on a much shorter timescale. This 
separation of timescales allows for an adiabatic elimina- 
tion of the variable r(x), assuming it to always be at the 
equilibrium point of Equation (27): 

r(x) = -L- 1 [Q r [z,z]]. (28) 

We remark that as ^t/lW < ^ f° r a U finite wave 
numbers k >0, the operator L r [r] is indeed invertible 
when excluding global translations in the set of possible 
perturbations of the trivial fixed point. From Equation 
(28), the coupled dynamics of OPM and RM is thus 
reduced to a third-order effective dynamics of the OPM: 

d t z(x) « L z [z] - (TIL; 1 [Q r [z, z]] , z] + N z 3 [z, z, z] 

N r 3 [z,z,z] (29) 

The nonlinearity N 3 [z,z,z] accounts for the coupling 
between OPM and RM. Its explicit analytical calculation 
for the EN model is rather involved and yields a sum 

12 

N r 3 [z,z,z] = ^4M[z,z,z]. 

i=i 

The individual nonlinear operators ^ are nonlinear 
convolution-type operators and are presented in the 
'Methods' section together with a detailed description of 
their derivation. Importantly, it turns out that the coeffi- 
cients j r are completely independent of the orientation 
stimulus ensemble. 

The adiabatic elimination of the retinotopic distortions 
results in an equation for the OPM (Equation (29)) 
which has the same structure as Equation (13), the only 
difference being an additional cubic nonlinearity. Due to 
this similarity, its stationary solutions can be determined 
by the same methods as presented for the case of a 
fixed retinotopy. Again, via weakly nonlinear analysis we 
obtain amplitude equations of the form Equation (15). 
The nonlinear coefficients gij and are determined 
from the angle-dependent interaction functions g(a) and 
f(oc). For the operator N r 3 [z,z,z], these functions are 
given by 

S» = 2(J 4 - G 2 e -2kW{cos «-!)) 

/r(«) = \ (gr{d) + gr{0C + 7t)) , 



verifying that, N 3 [z, z, z] is independent of the orienta- 
tion stimulus ensemble. Besides the interaction range ol 
A the continuity parameter r\ r e [0, oo] for the RM 
appears as an additional parameter in the angle-depen- 
dent interaction function. Hence, the phase diagram of 
the EN model will acquire one additional dimension 
when retinotopic distortions are taken into account. We 
note, that in the limit r\ r oo, the functions g r (a) and/ r 
(a) tend to zero and as expected one recovers the 
results presented above for fixed uniform retinotopy. 
The functions g r (oc) and f r (oc) are depicted in Figure 8 
for various interaction ranges <7/A and retinotopic conti- 
nuity parameters r\ r . 
Coupled essentially complex n-planforms 
In the previous section, we found that by an adiabatic 
elimination of the retinotopic distortions in the 
dynamics equations (26, 27) the system of partial inte- 
gro-differential equations can be reduced to a single 
equation for the OPM. In this case, the stationary solu- 
tions of the OPM dynamics are again planforms com- 
posed of a discrete set of Fourier modes 

N 

z(x) = J2 A i eikjX > ( 3 °) 

i 

with |k| = k c . However, each of these stationary plan- 
form OPM solutions induces a specific pattern of reti- 
notopic distortions by Equation (28). The joint mapping 
{X + r(x), z(x)} is then an approximate stationary solu- 
tion of Equations (26, 27) and will be termed coupled 
planform solution in the following. In contrast to other 
models for the joint mapping of orientation and visual 
space (e.g., [31,33,92]), the coupling between the repre- 
sentation of visual space and orientation in the EN 
model is not induced by model symmetries but a mere 
consequence of the joint optimization of OPM and RM 
that requires them to be matched to one another. 

For planforms given by Equation (30), it is possible to 
analytically evaluate Equation (28) and compute the 
associated retinotopic distortions r(x). After a somewhat 
lengthy calculation (see 'Methods' section), one obtains 

) ) (si) 

*(»(AjAjk) cos( A jfe x) + W{AjA k ) sin(Ajfcx)), 

with A Jk = kj - k /c and X r L (k) = -k 2 {r] r + e^V) . 
These retinotopic distortions represent superpositions of 
longitudinal modes (see Figure 3b). Hence, coupled 
planform stationary solutions of the EN dynamics do 
not contain any transversal mode components. Accord- 
ing to Equation (31), the pinwheel-free coupled 1-ECP 
state has the functional form {r(x) = 0, z(x) = A 0 e lklc }. 
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Figure 8 Angle-dependent interaction functions for the coupling between OPM and RM in the EN model, (a, b) g r {a) and f r {a) for r\ r = 
0.005 and a/A = 0.3 (a) and 0.1 (b). (c, d) g r {a) and f r {a) for r\ r = 0.05 and a/A = 0.3 (c) and 0.1 (d). (e, f) g r {a) and f r {a) for r\ r = 0.5 and a/A = 
0.3 (e) and 0.1 (f). 



This means that the OS solution does not induce any 
deviations from the perfect retinotopy as shown pre- 
viously from symmetry. This is not the case for the 
square pinwheel crystal (sPWC) 

£spwc(x) oc sm(k c Xi) + isin(fc c X2), 

the second important solution for undistorted retino- 
topy. Inserting this ansatz into Equation (31) and 

neglecting terms of order O ^£ _fe ^ 2 ^ \ Q r higher, we 

obtain 



e~ k ' a2 /fe c sin(2fe c xi)\ 
r5PWc(x)a a^(2^U c sin(2te)J- 

These retinotopic distortions are a superposition of 
one longitudinal mode in ^-direction and one in y- 
direction, both with doubled wave number ~ 2/c c . The 
doubled wave number implies that the form of retino- 
topic distortions is independent of the topological 
charge of the pinwheels. Importantly, the gradient of 
the retinotopic mapping R(x) = X + r sPWC (x) is 
reduced at all pinwheel locations. The coupled sPWC 
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is therefore in two ways a high coverage mapping as 
expected. First, the representations of cardinal and 
oblique stimuli (real and imaginary part of z(x)) are 
orthogonal to each other. Second, the regions of high- 
est gradient in the orientation map correspond to low 
gradient regions in the RM. 

In Figure 9, the family of coupled h-ECPs is displayed, 
showing simultaneously the distortions of the RM and 
the OPM. Retinotopic distortions are generally weaker 
for anisotropic h-ECPs and stronger for isotropic n- 
ECPs. However, for all stationary solutions the regions 
of high gradient in the orientation map coincide with 
low gradient regions (the folds of the grid) in the RM. 
This is precisely what is generally expected from a 
dimension-reducing mapping [21,62,63,91]. In the fol- 
lowing section, we will investigate which of these solu- 
tions become optimal depending on the two parameters 
olA and r\ r that parameterize the model. 

The impact of retinotopic distortions 

According to our analysis, at criticality, the nontrivial 
stable fixed points of the EN dynamics are determined 
by the continuity parameter 77 e (0, 1) for the OPM or, 
equivalently, the ratio a/ A = y/\og{l / rj) and the con- 
tinuity parameter r\ r for the mapping of visual space. 
We first tested for the stability of pinwheel-free orienta- 
tion stripe (OS) solutions and rPWC solutions of Equa- 
tion (15), with coupling matrices gy and/Jy as obtained 
from the nonlinearities in Equation (29). The angle 
which minimizes the energy [/ rPWC (Equation (20)) is 
not affected by the coupling between retinotopic and 
OPM and is thus again a = 7z74. By numerical evalua- 
tion of the criteria for intrinsic and extrinsic stability, 
we found both, OSs and sPWCs, to be intrinsically and 
extrinsically stable for all ol A and r\ r . 

Next, we tested for the stability of coupled n-ECP 
solutions for 2 < n < 20. We found all coupled w-ECP 
configurations with n > 2 to be intrinsically unstable for 
all ol A and r\ r . Evaluating the energy assigned to 
sPWCs and OSs, we identified two different regimes: (i) 
for shorter interaction range olA the sPWC is the mini- 
mal energy state and (ii) for larger interaction range ol 
A the optimum is an OS pattern as indicated by the 
phase diagram in Figure 10a. The retinotopic continuity 
parameter has little influence on the energy of the two 
fixed points. The phase border separating stripes from 
rhombs runs almost parallel to the 77 r -axis. We numeri- 
cally confirmed these analytical predictions by extensive 
simulations of Equation (3, 4) (see 'Methods' section for 
details). Figure 10c shows snapshots of a representative 
simulation with small interaction range (r = 0.1, olA = 
0.1 (77 = 0.67), r\ r = 77). After the initial symmetry break- 
ing phase, the OPM layout rapidly converges toward a 



crystalline array of pinwheels, the predicted optimum in 
this parameter regime (Figure 10c). Retinotopic devia- 
tions are barely visible. Figure 10b displays pinwheel 
density time courses for four such simulations. Note 
that in one simulation, the pinwheel density drops to 
almost zero. In this simulation, the OP pattern con- 
verges to a stripe-like layout. This is in line with the 
finding of bistability of rhombs and stripes in all para- 
meter regimes. Although the sPWC represents the glo- 
bal minimum in the simulated parameter regime, OSs 
are also a stable fixed point and, depending on the 
initial conditions, may arise as the final state of a frac- 
tion of the simulations. In the two simulations with pin- 
wheel densities around 3.4, patterns at later simulation 
stages consist of different domains of rhombic pinwheel 
lattices with a < 7z72. 

Figure 10d,e shows the corresponding analysis with 
parameters for larger interaction range r = 0.1, olA = 
0.15 (77 = 0.41), 77 r = 77. Here after initial pinwheel crea- 
tion, pinwheels typically annihilate pairwisely and the 
OPM converges to an essentially pinwheel-free stripe 
pattern, the predicted optimal solution in this parameter 
regime (Figure lOe). Retinotopic deviations are slightly 
larger. The behavior of the EN model for the joint opti- 
mization of RM and OPM thus appears very similar 
compared to the fixed retinotopy case. Perhaps surpris- 
ingly, the coupling of both feature maps has little effect 
on the stability properties of the fixed points and the 
resulting optimal solutions. 

As in the previous case, the structure of the phase dia- 
gram in Figure 10a appears somewhat counterintuitive. 
A high coverage and pinwheel-rich solution is the opti- 
mum in a regime with large OPM continuity parameter 
where discontinuities in the OPM such as pinwheels 
should be strongly penalized. A pinwheel-free solution 
with low coverage and high continuity is the optimum 
in a regime with small continuity parameter. As 
explained above, a large OPM continuity parameter at 
pattern formation threshold implies a small interaction 
range olA (see Equation (24)). In such a regime, the 
gain in coverage by representing many orientation sti- 
muli in a small area spanning the typical interaction 
range, e.g., with a pinwheel, is very high. Apparently this 
gain in coverage by a regular positioning of pinwheels 
outweighs the accompanied loss in continuity for very 
large OPM continuity parameters. This counterintuitive 
interplay between coverage and continuity thus seems to 
be almost independent of the choice of retinotopic con- 
tinuity parameters. 

The circular orientation stimulus ensemble contains 
only stimuli with a fixed and finite 'orientation energy' 
or elongation \s z \. This raises the question of whether 
the simple nature of the circular stimulus ensemble 
might restrain the dynamics of the EN model. The EN 
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Figure 9 Coupled n-ECPs as dimension-reducing solutions of the EN model Coupled n-ECP are displayed in visual space showing 
simultaneously the distortion of the RM and the OPM {a/A = 0.3 (77 = 0.028), r\ r = 77, circular stimulus ensemble). The distorted grid represents a 
the cortical square array of cells. Each grid intersection is at the receptive field center of the corresponding cell. Preferred stimulus orientations 
are color-coded as in Figure 2a. As in Figure 4, n and / enumerate the number of nonzero wave vectors and nonequivalent configurations of 
wave vectors with the same n, respectively. The coupled 1-ECP is a pinwheel-free stripe pattern without retinotopic distortion. Only the most 
anisotropic and the most isotropic coupled n-ECPs are shown for each n. Note that for all ECPs, high gradients within the orientation mapping 
coincide with low gradients of the retinotopic mapping and vice versa. Retinotopic distortions are displayed on a fivefold magnified scale for 
visualization purposes. 
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Figure 10 Phase diagram of the EN model with variable retinotopy for a circular stimulus ensemble [57,64-66,91]. (a) Regions of the r\ r - 
cr/A-plane in which n-ECPs or rPWCs have minimal energy, (b) Pinwheel density time courses for four different simulations of Equations (3, 4) 
with r = 0.1, a/ A = 0.13 (77 = 0.51), rj r = 77 (grey traces, individual realizations; red trace, mean value; black trace, realization shown in c). (c) OPMs 
(upper row), their power spectra (middle row), and RMs (lower row) obtained in a simulation of Equations (3, 4); parameters as in b. (d) Pinwheel 
density time courses for four different simulations of Equations (3, 4) with r = 0.1, cr/A = 0.3 (77 = 0.03), r\ r = 77 (grey traces, individual realizations; 
red trace, mean value; black trace, realization shown in e). (e) OPMs (upper row), their power spectra (middle row), and RMs (lower row) in a 
simulation of Equations (3, 4); parameters as in d. 



dynamics are expected to depend on the characteristics 
of the activity patterns evoked by the stimuli and these 
will be more diverse and complex with ensembles con- 
taining a greater diversity of stimuli. Therefore, we 
repeated the above analysis of the EN model for a richer 
stimulus ensemble where orientation stimuli are uni- 
formly distributed on the disk {s z , \s z \ < 2}, a choice 
adopted by a subset of previous studies, e.g., [19,25,81]. 
In particular, this ensemble contains unoriented stimuli 
with \s z \ = 0. Intuitively, the presence of these unor- 
iented stimuli might be expected to change the role of 
pinwheels in the optimal OPM layout. Pinwheels' popu- 
lation activity is untuned for orientation. Pinwheel cen- 
ters may therefore acquire a key role for the 
representation of unoriented stimuli. Nevertheless, we 
found the behavior of the EN model when considered 
with this richer stimulus ensemble to be virtually indis- 
tinguishable from the circular stimulus ensemble. Details 



of the derivations, phase diagrams and numerically 
obtained solutions are given in Appendix 1. 

Are there stimulus ensembles for which realistic, 
aperiodic maps are optimal? 

So far, we have presented a comprehensive analysis of 
optimal dimension-reducing mappings of the EN model 
for two widely used orientation stimulus distributions 
(previous sections and Appendix 1). In both cases, 
optima were either regular crystalline pinwheel lattices 
or pinwheel-free orientation stripes. These results might 
indicate that the EN model for the joint optimization of 
OPM and RM is per se incapable of reproducing the 
structure of OPMs as found in the visual cortex. Draw- 
ing such a conclusion is suggested in view of the appar- 
ent insensitivity of the model's optima to the choice of 
stimulus ensemble. The two stimulus ensembles consid- 
ered so far however do not exhaust the infinite space of 
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stimulus distributions that are admissible in principle. 
From the viewpoint of 'biological plausibility' it is cer- 
tainly not obvious that one should strive to examine sti- 
mulus distributions very different from these, as long as 
the guiding hypothesis is that the functional architecture 
of the primary visual cortex optimizes the joint repre- 
sentation of the classical elementary stimulus features. 
If, however, stimulus ensembles were to exist, for which 
optimal EN mappings truly resemble the biological 
architecture, their characteristics may reveal essential 
ingredients of alternative optimization models for visual 
cortical architecture. 

Adopting this perspective raises the technical question 
of whether an unbiased search of the infinite space of 
stimulus ensembles only constrained by the model's 
symmetries (Equations (5) to (9)) is possible. To answer 
this question, we examined whether the amplitude equa- 
tions (15) can be obtained for an arbitrary orientation 
stimulus distribution. Fortunately, we found that the 
coefficients of the amplitude equations are completely 
determined by the finite set of moments of order less 
than 5 of the distributions. The approach developed so 
far can thus be used to comprehensively examine the 
nature of EN optima resulting for any stimulus distribu- 
tion with finite fourth-order moment. While such a 
study does not completely exhaust the infinite space of 
all eligible distributions, it appears to only exclude 
ensembles with really exceptional properties. These are 
probability distributions with diverging fourth moment, 
i.e., ensembles that exhibit a heavy tail of essentially 
'infinite' orientation energy stimuli. 

Since the coupling between OPM and RM did not 
have a large impact in the case of the two classical sti- 
mulus ensembles, we start the search through the 
space of orientation stimulus ensembles by considering 
the EN model with fixed retinotopy r(x) = 0. The coef- 
ficients a t for the nonlinear operators Nf[z,z,z\ in 
Equation (14) for arbitrary stimulus ensembles are 
given by 

"! 16a 6 2a 4 + la 1 " 2 8;rff 6 32jra 8 "3 64jT(T 8 + 16na 6 

, _ (ki 4 ) . (Wi 2 ) i , (i^i 4 ) , , , 

n ISzl 4 (\SA 2 ) n |5,| 4 3 |S Z [ 4 W^) 

dl = 48^10 - fl 8 = g^vTo ^9 = " 256^%" 

l* Z | 4 (ISzl 2 „ Hzl 4 

The corresponding angle-dependent interaction func- 
tions are given by (see 'Methods' section) 

m - M (l - 2e~^ - (l - 2e -^«-)) 

(i s pi , " 6 iA (33) 

f(a) = ^1 (l - e- 2k ' a (cosh(2fe2 CT 2 cosa) + 2 cosh (fe2 CT 2 cos a)) + 2e~% a ) 
+ ^ (e- 2k > 2 cosh(2fe2 CT 2 cosa) _ ^ + fcp e -2fe^ 2 sinh 4 ( 1/2 fe2 ff 2 CQSa ) 



Again, without loss of generality, we set (|s z | 2 ) = 2. At 
criticality, both functions are parameterized by the conti- 
nuity parameter 77 e (0, 1) for the OPM or, equivalently, 
the interaction range a/ A = ^yiogfTT^) and the fourth 
moment (|s z | 4 ) of the orientation stimulus ensemble. The 
fourth moment, is a measure of the peakedness of a sti- 
mulus distribution. High values generally indicate a 
strongly peaked distribution with a large fraction of non- 
oriented stimuli (\s z \ 4 « 0), together with a large fraction 
of high orientation energy stimuli {\s z \ 4 large). 

The dependence of g(a) on the fourth moment of the 
orientation stimulus distribution and f(a) suggests that 
different stimulus distributions may indeed lead to dif- 
ferent optimal dimension-reducing mappings. 

The circular stimulus ensemble possesses the minimal 
possible fourth moment, with (|s z | 4 ) = «|s z | 2 )) 2 = 4. The 
fourth moment of the uniform stimulus ensemble is (| 
5 Z | 4 ) = 16/3. The angle-dependent interaction functions 
for both ensembles (Equation (25), Figure 22 in Appen- 
dix 1) are recovered, when inserting these values into 
Equation (33). 

To simplify notation in the following, we define 

54 = (^| 4 >-(l^| 2 ) 2 = (l^| 4 >- 4 

as the parameter characterizing an orientation stimu- 
lus distribution. This parameter ranges from zero for 
the circular stimulus ensemble to infinity for ensembles 
with diverging fourth moments. Figure 11 displays the 
angle-dependent interaction functions for different 
values of o7A and s 4 . In all parameter regimes, g(a) and 
f(a) are larger than zero. The amplitude dynamics are 
therefore guaranteed to converge to a stable stationary 
fixed point and the bifurcation from the nonselective 
fixed point in the EN model is predicted to be supercri- 
tical in general. 

By evaluating the energy assigned to the rPWC and n- 
ECPs, we investigated the structure of the two-dimen- 
sional phase space of the EN model with an arbitrary 
orientation stimulus distribution. First, it is not difficult 
to show that the angle a which minimizes the energy 
^rpwc (Equation (20)) of an rPWC is a = n/4 for all ol 
A and s 4 . Hence, a square lattice of pinwheels (sPWC) 
is in all parameter regimes energetically favored over 
any other rhombic lattice configuration of pinwheels. 
Figure 12 displays the phase diagram of the EN model 
with an arbitrary orientation stimulus distribution. For 
orientation stimulus distributions with small fourth 
moments, optimal mappings consist of either parallel 
pinwheel-free stripes or quadratic pinwheel crystals. 
These distributions include the circular and the uniform 
stimulus ensembles with s 4 = 0 and s 4 = 4/3. Above a 
certain value of the fourth moment around s 4 = 6, n- 
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Figure 1 1 Angle-dependent interaction functions for the EN model with fixed retinotopy for different fourth-moment values of the 

orientation stimulus distribution and effective interaction-widths, (a, b) g{a) and f{a) for s 4 = 8 and cr/A = 0.1 (a) and cr/A = 0.3 (b). (c, d) 

g(a) and f{a) for s 4 = 20 and cr/A = 0.1 (c) and cr/A = 0.3 (d). (e, f) g(a) and f{a) for s 4 = 100 and cr/A = 0.1 (e) and cr/A = 0.3 (f). 
I J 



ECPs with h >2 become optimal mappings. For a short 
interaction range cr/A, hexagonal pinwheel crystals dom- 
inate the phase diagram in a large region of parameter 
space. With increasing interaction range, we observe a 
sequence of phase transitions by which higher n-ECPs 
become optimal. For n >3, these optima are spatially 



aperiodic. In all parameter regimes, we found that the 
n-ECP with the most anisotropic mode configuration 
(Figure 4c, left column) is the energetically favored state 
for n >3. Pinwheel densities of these planforms are indi- 
cated in Figure 12 and are typically smaller than 2.0. 
We note that this is well below experimentally observed 
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Figure 12 Stripe-like, crystalline, and quasi-crystalline cortical representations as optimal solutions to the mapping of orientation 
preference with fixed uniform retinotopy in the EN model. The graph shows the regions of the s 4 -cr/A-plane in which n-ECPs or sPWCs 
have minimal energy. For n > 3, pinwheel densities of the energetically favored n-ECP configuration are indicated. 



pinwheel density values [38]. Optimal mappings of 
orientation preference for finite fourth moment in the 
EN model are thus either orientation stripes, periodic 
arrays of pinwheels (hexagonal, square), or aperiodic 
pinwheel arrangements with low pinwheel density. 



We numerically tested these analytical predictions by 
simulations of Equation (3) (r(x) = 0) with two addi- 
tional stimulus ensembles with s 4 = 6 and s 4 = 8 (see 
'Methods' section). Figure 13a shows snapshots of a 
simulation with (r = 0.1, cr/A = 0.2 (77 = 0.2)) and s 4 = 6 
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Figure 13 Approaching crystalline n-ECP optima in the EN model with fixed renotopy. (a) OPMs (upper row) and their power spectra 
(lower row) in a simulation of Equation (3) with r(x) = 0, r = 0.1, o/A = 0.2 and s 4 = 6 (see also Additional file 3). The predicted optimum is the 
2-ECP (black frame), (b) Pinwheel density time courses for four different simulations (parameters as in a; gray traces, individual realizations; black 
trace, simulation in a; red trace, mean value), (c) OPMs (upper row) and their power spectra (lower row) in a simulation of Equation (3) with r(x) 
= 0, r = 0.1, o/A = 0.3 and s 4 = 8 (see also Additional file 4). The predicted optimum is the anisotropic 3-ECP (black frame), (d) Pinwheel density 
time courses for four different simulations (parameters as in c; gray traces, individual realizations; black trace, simulation in c; red trace, mean 
value). 



(see also Additional file 3). After the initial phase of pat- 
tern emergence, the OPM layout converges toward an 
arrangement of fractured stripes which resembles the 2- 
ECP state (Figure 13a, most right), the optimum pre- 
dicted in this regime. In the power spectra, two distinct 
peaks of the active modes are clearly visible in the final 
stages of the simulation (Figure 13a, lower row). The 2- 
ECP state is exotic in the sense that it is the only n-ECP 
containing line defects and thus the pinwheel density is 
not a well-defined quantity. This explains the pro- 
nounced numerical variability in the measured pinwheel 
densities in simulations during the convergence toward 
a 2-ECP state (Figure 13b). 

Figure 13c shows snapshots of a simulation with (r = 
0.1, ct/A = 0.2 (77 = 0.2)) and s 4 = 8, Gaussian stimulus 
ensemble) (see also Additional file 4). After the initial 
phase of pattern emergence, the OPM layout converges 
toward a regular hexagonal arrangement of pinwheels 
which resembles the anisotropic 3-ECP (Figure 13c, far 
right), the optimum predicted in this regime. In the 
power spectra, three distinct peaks forming an angle of 
60 degrees are clearly visible in the later stages of the 



simulation (Figure 13c, lower row). Pinwheel densities in 
the simulations consistently approach the theoretically 
predicted value of 2 cos(7r/6) = 1.73 (Figure 13d). 
Permutation symmetric limit 

In the previous section, we uncovered a parameter 
regime for the EN model in which optimal solutions are 
spatially aperiodic. This can be viewed as a first step 
toward realistic optimal solutions. In the identified 
regime, however, among the family of w-ECPs only 
those with pinwheel densities well below experimentally 
observed values [38] are energetically favored (see Figure 
12). In this respect, the repertoire of aperiodic optima of 
the EN model differs from previously considered 
abstract variational models for OPM development 
[35,36,38,39]. In these models, an energetic degeneracy 
of aperiodic states with low and high pinwheel densities 
has been found which leads to a pinwheel statistics of 
the repertoire of optimal solutions that quantitatively 
reproduces experimental observations [38,93]. What is 
the reason for this difference between the two models? 
In [35], the energetic degeneracy of aperiodic states with 
low and high pinwheel densities was derived from a so- 
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called permutation symmetry 

N z 3 [u, v, w] = N z 3 [w, u, v], (34) 

of the cubic nonlinearities of the model. It can be 
easily seen, that the cubic nonlinearities obtained in the 
third order expansion of the EN model do not exhibit 
this permutation symmetry (see 'Methods' section). As 
shown by Reichl [94], the absence of permutation sym- 
metry can lead to a selection of a subrange of pinwheel 
densities in the repertoire of optima of OPM models. 
Depending on the degree of permutation symmetry 
breaking, the family of optima of such models, albeit 
encompassing aperiodic OPM layouts, may consist of 
layouts with either unrealistically low or high pinwheel 
densities. Furthermore, for very strong permutation 
symmetry breaking, stationary solutions from solution 
classes other than the h-ECPs and rPWCs with low or 
high pinwheel densities may become optima of models 
for OPM development. In order to determine a regime 
in which the EN model optima quantitatively resemble 
experimentally observed OPM layouts, it is therefore 
important to quantify the degree of permutation sym- 
metry breaking in the EN model and to examine 
whether permutation symmetric limits exist. As shown 
in the 'Methods' section, any cubic nonlinearity 
N3 [z, z, z] that obeys Equation (34) has a corresponding 
angle-dependent interaction function g(a) which is n- 
periodic. Therefore, we examine the degree of permuta- 
tion symmetry breaking in the EN model by comparing 
the angle-dependent interaction function g(a) of its 
third order expansion (see Equation 33 and Figure 11) 
to the 7r-periodic function g pm (oc) = 1/2 (g(cc) + g(a + 



n)). This 'permutation-symmetrized' part of the angle- 
dependent interaction function of the EN model for 
general orientation stimulus ensembles reads 

gpm{a) = l^l e - 2k ^ 2 smh 4 {l/2k 2 ycosa) 

~ ^pT e ~ 2kW ((cosh(2^a 2 cosa) - 2 cosh a 2 cos a)) -2^ - e 2 ^) (35) 
+ 2^2 (l + e ~ lKal cosh(2fe 2 a 2 cos a)) . 

A comparison between g pm {a) and g(a) is depicted in 
Figure 14a-d. It shows that essentially insensitive to the 
interaction range a/A, at large values of the fourth 
moment original and permutation symmetrized angle- 
dependent interaction functions converge. We quanti- 
fied the degree of permutation symmetry breaking with 
the parameter 

d= ll^-^Jl2 

II gh 

This parameter is zero in the case of a permutation 
symmetric cubic nonlinearity. In the case of a g-function 
completely antisymmetric around a = n/2, the para- 
meter is either plus or minus one, depending on 
whether the maximum of g pm is at zero or n. If d is 
smaller than zero, low pinwheel densities are expected 
to be energetically favored and vice versa. The values of 
d in parameter space is depicted in Figure 14e. It is 
smaller than zero in the entire phase space, implying a 
tendency for low pinwheel density optimal states, in 
agreement with the phase diagram in Figure 12. Permu- 
tation symmetry breaking is largest for a/A around 0.25 
and small fourth moment values of the orientation sti- 
mulus distribution. It decays to zero for large fourth 
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Figure 14 Quantifying permutation symmetry breaking in the EN model, (a-d) g{a) (red traces) and the 'permutation symmetrized' function 
g pm {a) = ]/2{g{a) + g{a + n)) (blue traces, see Equation (35)) for a/ A = 0.1 and 0.3 and s 4 = 6 and 100. (e) Permutation symmetry parameter d 
(Equation 36) in the EN model with fixed retinotopy. Permutation symmetry breaking is largest for cr/A ~ 0.25 and small s 4 . In the limit s 4 °° 
permutation symmetry is restored. 
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moments proportionally to l/s 4 as can be seen by insert- 
ing Equations (33) and (35) into Equation (36). In the 
infinite fourth moment limit s 4 — > the cubic nonli- 
nearities of the third-order expansion of the EN model 
become permutation symmetric. 

In this case, the EN model is parameterized by only 
one parameter, the effective intracortical interaction 
range a/A and we obtain a rather simple phase diagram 
(Figure 15). Optimal solutions are n-ECPs for increasing 
a/ A and we observe a sequence of phase transitions 
toward a higher number of active modes and therefore 
more complex spatially aperiodic OPM layouts. Impor- 
tantly, for a subregion in the phase diagram with given 
number of active modes, all possible n-ECP mode con- 
figurations are energetically degenerate. It is precisely 
this degeneracy that has been previously shown to result 
in a pinwheel statistics of the repertoire of aperiodic 
optima which quantitatively agrees with experimental 
observations [38]. Therefore, our unbiased search in fact 
identified a regime, namely a very large effective interac- 
tion range and infinite fourth moment of the orientation 
stimulus ensemble, in which the EN model formally pre- 
dicts which quantitatively reproduce the experimentally 
observed VI architecture. 

Unexpectedly, however, this regime coincides with the 
limit of applicability of our approach. Permutation sym- 
metry is exactly obtained by approaching stimulus distri- 
bution with diverging fourth moment for which the 
amplitude equations may become meaningless. We 
would generally expect that the EN for very large but 
finite fourth moment can closely resemble a permuta- 
tion symmetric model. However, to consolidate the 



relevance of this regime, it appears crucial to establish 
the robustness of the limiting behavior to inclusion of 
retinotopic distortions. 

Optimal solutions of the EN model with variable 
retinotopy and arbitrary orientation stimulus ensembles 

In the EN model for the joint mapping of visual space 
and orientation preferences, the angle-dependent inter- 
action functions depend on four parameters: 77, a, the 
fourth moment (|s z | 4 )of the stimulus ensemble and r] r . 
By setting a = a* (77), we are left with three free para- 
meters at criticality. Therefore, a three-dimensional 
phase diagram now completely describes pattern selec- 
tion in the EN model. For better visualization, in Figure 
16, we show representative cross sections through this 
three-dimensional parameter space for fixed r\ r . First, we 
note the strong similarity between the phase diagram 
for fixed retinotopy (Figure 12) and the cross sections 
through the phase diagrams for the joint mappings 
shown in Figure 16. This expresses the fact that retino- 
topic mapping and OPM are only weakly coupled or 
mathematically, g r (a) « g(a) in all parameter regimes 
(see Appendix 2). Again, for distributions with small 
fourth moment, optimal mappings consist of either pin- 
wheel-free orientation stripes or sPWCs. Above a cer- 
tain fourth moment value around s 4 = 6, higher coupled 
n-ECPs are optimal. For small interaction range o7A, 
hexagonal pinwheel crystals (coupled 3-ECPs) represent 
optimal mappings in a large fraction of parameter space. 
With increasing o7A, we observe a sequence of phase 
transitions by which higher n-ECPs become optimal. 
Anisotropic planforms at the lower end of the spectrum 
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Figure 15 Phase diagram of the EN model with fixed retinotopy in the permutation symmetric limit s 4 -> 00. The graphs show the 
regions on the o/A-axis (lower axis) and the corresponding 77-axis (upper axis), where n-ECPs or sPWCs have minimal energy. High n-ECPs (n £ 
10) exhibit universal pinwheel statistics. Note however the extremely small 77-values for large a/A 
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of pinwheel densities are always energetically favored 
over high pinwheel density layouts. The only difference 
between the cross sections is that the region covered by 
sPWCs increases for decreasing r\ r . The phase diagram 
for large r\ r - 10 is virtually indistinguishable from the 
phase diagram in Figure 16. 

Optimal mappings of orientation preference are thus 
either orientation stripes, periodic arrays of pinwheels 
(hexagonal, quadratic) or quasi-periodic pinwheel arrays 
with low pinwheel density. Retinotopic distortions lead 
to lower gradients of the retinotopic mapping at high 
gradient regions of the OPM. This is in line with some 
of the experimental evidence [55,95] but contradicts 
others [96]. 

Most importantly, we note that the results on permu- 
tation symmetry breaking in the fixed retinotopy case 
are not altered by allowing for retinotopic distortions. 
Since g r (cc) does not depend on the fourth moment of 
the orientation stimulus distribution, non-permutation 
symmetric terms decay as l/s 4 for large s 4 . 

Hence, in the limit s 4 — » °o, permutation symmetry is 
restored and we recover the phase diagram in Figure 
15 also for the EN model with variable retinotopy 
independent of 7] r . As the energy contribution of reti- 
notopic deviations r(x) becomes negligible in the infi- 
nite fourth moment limit, the optima are then simply 
the corresponding coupled n-ECPs and these states are 
energetically degenerate for fixed n. For very large 
effective interaction range and infinite fourth moment 
of the orientation stimulus ensemble, the EN model 
with variable retinotopy is able to quantitatively repro- 
duce the experimentally observed pinwheel statistics in 
OPMs. It furthermore predicts reduced gradients of 
the visual space mapping at high gradient regions of 
the OPM. 

Finite stimulus samples and discrete stimulus ensembles 

Our reexamination of the EN model for the joint opti- 
mization of position and orientation selectivity has been 
so far carried out without addressing the apparently fun- 
damental discrepancy between our results and the large 
majority of previous reports. Since the seminal publica- 
tion of Durbin and Mitchison [21], numerous studies 
have used the EN model to simulate the development of 
visual cortical maps or to examine the structure of opti- 
mal mappings by numerical simulation 
[58,62-65,79,97,98]. These studies have either used the 
circular or the uniform orientation stimulus ensemble 
for which, to the best of our knowledge, the only two 
nontrivial stationary solutions are square pinwheel crys- 
tals or orientation stripes. Furthermore, we found that 
the gradient descent dynamics seems to readily converge 
to the respective minima of the EN free energy. This 
indicates that other local minima and more complex 



intrinsically aperiodic states are not dominant in this 
model. In fact, we found that all aperiodic stationary 
solutions we could perturbative calculate analytically are 
unstable and thus represent hyperbolic saddle points 
and not local minima. As these stable solutions barely 
resemble experimentally observed OPMs, it is not 
obvious how the EN model in all of these studies could 
appear as a model well suited to describe the complex 
layout of real cortical orientation maps. Prior studies 
however often used computational methods different 
from our fixed parameter steepest descent simulations. 

Two alternative approaches have been used predomi- 
nantly to study dimension reducing mappings for corti- 
cal representations. These methods have been applied to 
both the EN model and the other widely used dimen- 
sion reduction model, the self-organizing feature map 
(SOFM), originally introduced by Kohonen [59]. The 
simplest way to compute mappings from a high dimen- 
sional feature space onto the two-dimensional model 
cortex is by iterating the following procedure for a large 
number of randomly chosen stimuli (e.g., 
[56,57,66-68,99,100]): (i) Stimuli are chosen one at a 
time randomly from the complete feature space, (ii) The 
activation function for a particular stimulus is com- 
puted. In the case of the EN model, this activation func- 
tion can acquire a rather complex form with multiple 
peaks (see 'Discussion' section). In the case of an 
SOFM, this activation function is a 2D-Gaussian. (iii) 
The preferred features of the cortical grid points are 
updated according to a discretized version of Equations 
(3, 4) or the corresponding equations for the SOFM 
model. Typically, this procedure is repeated on the 
order of 10 6 times. The resulting layout is then assumed 
to at least approximately solve the dimension reduction 
problem. In many studies, small stimulus sets have been 
chosen presumably for computational efficiency and not 
assuming specifically that the cortex is optimized for a 
discrete finite set of stimuli. In [21] for instance, a set of 
216 stimuli was used, that was likely already at the limit 
of computing power available at this time. 

In a more refined approach, the EN model as well as 
Kohonen's SOFM model have been trained with a finite 
set of stimuli (typically with on the order of 10 3 to 10 4 ) 
and the final layout of the model map has been obtained 
by deterministic annealing [101], i.e., by gradually redu- 
cing the numerical value of o in a numerical minimiza- 
tion procedure for the energy functional T at each 
value of <7 (see e.g., [21,64,65,79] and see 'Methods' sec- 
tion). In such simulations, often nonperiodic boundary 
conditions were used. One might suspect in particular 
the second approach to converge to OPM layouts 
deviating from our results. It is conceivable in principle, 
that deterministic annealing might track stationary solu- 
tions across parameter space that are systematically 
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Figure 16 Stripe-like, crystalline, and quasi-crystalline cortical representations as optimal solutions to the joint mapping problem of 
visual space and orientation preference in the EN. (a-d) Phase diagrams for the joint mapping of visual space and orientation preference in 
the EN near criticality for rj r = 0 (a), r\ r = 0.01 (b), r\ r = 0.1 (c), and r\ r = 10 (d). The graphs show the regions of the s 4 -o/A-plane in which 
coupled n-ECPs or sPWCs have minimal energy. For n > 3, pinwheel densities of the energetically favored n-ECP configuration are indicated. 
Note the strong similarity between the phase diagrams and the phase diagrams in the fixed retinotopy case (Figure 12). 



missed by both, our continuum limit analytical calcula- 
tions as well as our descent numerical simulations. 

To assess the potential biases of the different 
approaches, we implemented (i) finite stimulus sampling 



in our gradient descent simulations and (ii) studied the 
results of deterministic annealing simulations varying 
both the size of the stimulus set as well as the type of 
boundary conditions applied. 
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We simulated Equations (3, 4) with finite sets of sti- 
muli of different sizes (see 'Methods' section), drawn 
from the circular stimulus ensemble. Following [21,25], 
7] was set to a small value (77 = 0.025) such that the 
optimal configuration for the joint mapping of visual 
space and orientation preference is the coupled 1-ECP 
(see Figure 10), i.e., a pattern of parallel orientation 
stripes without any retinotopic distortion (see Figure 9). 
Figure 17 displays representative simulations for stimu- 
lus sets of size N = 216 (as used in [21]) (a), N = 10 5 
(b), N = 10 6 (c) stimuli. Simulation time t is measured 
in units of the intrinsic time scale r (see 'Methods' sec- 
tion). For N = 216 stimuli, RM and OPM quickly reach 
an apparently stationary configuration with a large num- 
ber of pinwheels at around t = 20r. Power is distributed 
roughly isotropically around the origin of Fourier space 



(k = 0). The stable OPM lacks a typical length scale 
and, expressing the same fact, the power spectrum lacks 
the characteristic ring of enhanced Fourier amplitude. 
Retinotopic distortions are fairly pronounced. Both 
obtained maps resemble the configurations reported in 
[21]. 

For N = 10 5 stimuli, we find that OPMs exhibit a 
characteristic scale (see dark shaded ring in the power 
spectrum) and a dynamic rearrangement of the maps 
persists at least until t = 200r. Stripe-like OP domains 
are rapidly generated via pairwise pinwheel annihilation 
for t >10r. Retinotopic distortions are fairly weak. For N 
= 10 6 stimuli, again OPMs exhibit a characteristic scale 
(see dark shaded ring in the power spectrum) and the 
map dynamics persists beyond t = 200r. A larger frac- 
tion of the pinwheels annihilate pairwisely compared to 
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Figure 17 Development of OPM and retinotopic distortions in EN simulations with fixed stimulus sets of different sizes, (a) OPMs (left), 
their power spectra (middle) and RMs (right) for f = 1 0r (upper row) and t = 300r (lower row) obtained in simulations with fixed stimulus set (77 
= 0.028, a/A = 0.3, s4 = 4/3, 216 stimuli), (b) 10 5 stimuli (all other parameters as in a), (c) 10 6 stimuli (all other parameters as in a). Large stripe- 
like OP domains are generated via pairwise pinwheel annihilation for large simulation times. Retinotopic distortions are fairly weak.(d) Pinwheel 
density time course for EN simulations with fixed stimulus sets of different sizes, including the simulations from a to c (red, green, blue traces 
216, 10 5 , 10 6 stimuli) (all other parameters as in a). Dashed lines represent individual simulations, solid lines an average over four simulations. 
Note, that the pinwheel density rapidly decays below 2.0 in both cases, and in particular for 10 6 stimuli, the OPM pattern acquires large stripe- 
like regions. 
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N = 10 5 stimuli, leading progressively to a pattern with 
large stripe-like domains. Retinotopic distortions are 
fairly weak. For both cases with massive stimulus sam- 
pling (N = 10 5 , N = 10 6 ), the pinwheel density rapidly 
drops below the range observed in tree shrews, galagos 
and ferrets and than further decreases during subse- 
quent map rearrangement. In summary, the more sti- 
muli are chosen for the optimization procedure, the less 
pinwheels are preserved in the pattern of orientation 
preference and the more the resulting map resembles 
the analytically obtained optimal solution. Deterministic 
annealing approaches which change parameters of the 
energy functional during the computational minimiza- 
tion process differ more fundamentally from our gradi- 
ent descent simulations than the iterative schemes used 
with fixed parameters. Studies using deterministic 
annealing in addition frequently used nonperiodic 
boundary conditions (e.g. [64,65,79]). To study all 
potential sources of deviating results, we implemented 
deterministic annealing for the EN energy function (see 
'Methods' section, Equation (46)) for periodic bound- 
aries, nonperiodic boundary conditions as well as ran- 
dom and grid-like finite stimulus ensembles (see 
'Methods' section). We closely follow the refined meth- 
ods used in [64,65,79] and performed deterministic 
annealing simulations for the EN model with retinotopic 
distortions and stimuli drawn from the circular stimulus 
ensemble. 

Figures 18a and 19a display representative simulations 
for random stimulus sets of size N = 10 3 , N = 10 4 and 
N = 10 5 for periodic boundary conditions (Figures 18a) 
and nonperiodic boundary conditions (Figures 19a). 
Furthermore depicted are the pinwheel densities of sta- 
tionary solutions as well as their energies, relative to the 
energy of a pinwheel-free stripe solution (see 'Methods' 
section) for different annealing rates <f (Figures 18b-d 
and 19b-d). Figures 18e-g and 19e-g additionally show 
the statistics of nearest neighbor (NN) pinwheel dis- 
tances as well as the SD of the pinwheel densities for 
randomly selected subregions in the OPM as introduced 
in [38], averaged over four simulations with N = 10 5 . To 
facilitate comparison, we superimposed fits to the 
experimentally observed statistics [38] for orientation 
maps in tree shrews, ferrets and galagos. 

When annealing with periodic boundary conditions, 
the maps found with deterministic annealing essen- 
tially resemble our gradient descent dynamics simula- 
tions. The larger the set of stimuli, the more stripe-like 
are the OPMs obtained (Figure 18a,b). Furthermore, 
the more carefully we annealed, the lower the pinwheel 
density of the obtained layouts (Figure 18c). For N = 
10 5 , the pinwheel density averaged over four simula- 
tions with annealing rate 0.999 was p = 2.04 As 
expected, the energy of the final layouts decreased 



with slower annealing rates (Figure 18d). However, 
when starting from random initial conditions, the 
energy of the final layouts found was always higher 
compared to the energy of a pinwheel-free stripe solu- 
tion (see 'Methods' section for details), which is the 
predicted optimum for the circular stimulus ensemble. 
NN-pinwheel distance histograms are concentrated 
around half the typical column spacing and in particu- 
lar pinwheel pairs with short distances are lacking 
completely (Figure 18e,f). For nonperiodic boundary 
conditions and random stimuli, we found that retino- 
topic distortions are more pronounced than for peri- 
odic boundary conditions. They however decreased 
with increasing number of stimuli. For large the stimu- 
lus numbers, we observed stripe-like orientation pre- 
ference domains which are interspersed with lattice- 
like pinwheel arrangements (see Figure 19c), lower 
row, upper left corner of the OPM). For N = 10 5 , the 
pinwheel density averaged over four simulations with 
annealing rate 0.999 was p = 2.71. 

Similarly to the results for periodic boundary condi- 
tions, short distance pinwheel pairs occur less frequently 
than in the experimentally observed maps, indicating an 
increased regularity in the pinwheel distances compared 
to real OPMs (Figure 19e,f). This regularity is further 
indicated by a smaller exponent of the SD compared to 
the Poisson process (Figure 19g). The perfect stripe-like 
solution is not the optimum for nonperiodic boundaries. 
The energy of the map layouts found with very slow 
annealing rates is slightly lower than the energy of the 
pinwheel-free OPM layout (Figure 19d). We note that 
the layout of the OPM at the boundaries does not differ 
substantially from the layout inside the simulated 
domain, suggesting that boundary effects affect the 
entire simulated domain for the relatively small region 
treated. Finally, we performed simulations with grid-like 
stimulus patterns as e.g., used in [64,65]. These simula- 
tions displayed a strong tendency toward rhombic pin- 
wheel arrangements, i.e., the second stable stationary 
solution found for the circular stimulus ensemble. We 
refer to Appendix 3 for further details. In summary, our 
results for the discrete EN model with deterministic 
annealing largely agree with the analytical results. Irre- 
spective of the numerical methodology, the emerging 
map structure for large numbers of stimuli is confined 
to the states predicted by our analytical treatment of the 
continuum formulation of the EN. This behavior is 
expected because the energies underlying the determi- 
nistic annealing and the steepest descent simulations are 
mathematically equivalent (see 'Methods' section). In 
any kind of deterministic annealing simulation we 
tested, resulting patterns were patchworks of the two 
fundamental stable solutions identified by the analytical 
treatment: pinwheel free stripes and square lattices of 
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Periodic boundary conditions 
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Figure 18 The EN model with periodic boundary conditions, solved with deterministic annealing, (a) OPMs (left) and RMs (right) for N = 
10 3 (upper row), N = 10 4 (middle row) and N = 10 5 (lower row) random stimuli and periodic boundary conditions (annealing rate % = 0.999, see 
'Methods' section). f3 is the continuity parameter in the conventional definition of the EN model (see 'Methods' section, Equation 46) and is 
scaled, such that a comparable number of columns is emerging in the simulations for each size of the stimulus set. (b) Pinwheel densities of EN 
solutions for different numbers of stimuli, % = 0.999. (c) Pinwheel densities of EN solutions for 10 5 stimuli and different annealing rates, (d) 
Energies of solutions for 10 5 stimuli, relative to the energy of a pinwheel-free stripe solution (see 'Methods' section) for different annealing rates, 
(b-d) Crosses mark individual simulations, red line indicates average values, (e, f) Statistics of nearest neighbor pinwheel distances for pinwheels 
of (e) arbitrary and (f) opposite and equal charge for 10 5 random stimuli and periodic boundary conditions, averaged over four simulations (red 
curves). Black curves represent fits to the experimental data from [38]. (g) Standard deviations (SD) of pinwheel densities estimated from 
randomly selected regions in the OPM. Black dashed curve indicates SD for a two-dimensional Poisson process of equal density. 



pinwheels. Such patchworks are spatially more compli- 
cated than perfect stripes or crystals. Nevertheless, they 
qualitatively differ in numerous respects from the 
experimentally observed spatial arrangements (see Fig- 
ures 18, 19, and 28 in Appendix 3). How the fundamen- 
tal stable solutions are stitched together somewhat 
differs between the different kinds of simulations. For 
instance, using a grid-like stimulus ensemble with non- 
periodic boundary conditions apparently energetically 
favors the rPWC compared to the pinwheel-free stripe 
regions (see Figure 27 in Appendix 3). In summary, 
while some of the patterns obtained by deterministic 



annealing might be called 'good-looking' maps, all of 
them substantially deviate from the characteristics of 
experimentally observed pinwheel arrangements. 

We conclude that the differences between our results 
and those of previous studies are most likely due to the 
small finite stimulus samples used largely for reasons of 
computational tractability. Deterministic annealing using 
stimulus samples that fill the feature space converges to 
the same types of patterns found by perturbation theory. 
We further conclude that our methods do not systema- 
tically miss biologically relevant local minima of the 
classical EN energy function. 
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Non-periodic boundary conditions 
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Figure 19 The EN model with nonperiodic boundary conditions, solved with deterministic annealing, (a-g) As Figure 18, but for 
nonperiodic boundary conditions. 



Discussion 

Summary 

In this study, we examined the solutions of what is per- 
haps the most prominent optimization model for the 
spatial layout of orientation and RMs in the primary 
visual cortex, the EN model We presented an analytical 
framework that enables us to derive closed-form expres- 
sions for hyperbolic fixed points, local and global 
minima, and to analyze their stability properties for arbi- 
trary optimization models for the spatial layout of 
OPMs and RMs. Using this framework, we systemati- 
cally reexamined previously used instantiations of the 
EN model, dissecting the impact of stimulus ensembles 
and of interactions between the two maps on optimal 
map layouts. To our surprise, the analysis yielded vir- 
tually identical results for all of these model instantia- 
tions that substantially deviate from previous numerical 
reports. Pinwheel-free orientation stripes and crystalline 



square lattices of pinwheels are the only optimal dimen- 
sion-reducing OPM layouts of the EN model. Both 
states are generally stable but exchange their roles as 
optima and local minima at a phase border. Numerical 
simulations of the EN gradient descent dynamics as well 
as simulations utilizing deterministic annealing con- 
firmed our analytical results. For both processes, the 
initially spatially irregular layouts rapidly decayed into a 
patchwork of stripe-like or crystal-like local regions that 
then became globally more coherent on longer time- 
scales. Pinwheel-free solutions were approached after an 
initial phase of pattern emergence by pairwise pinwheel 
annihilation. Crystalline configurations were reached by 
the generation of additional pinwheels and pinwheel 
annihilation together with a coordinated rearrangement 
toward a square lattice. These results indicate that lay- 
outs which represent an optimal compromise of cover- 
age and continuity for retinotopy and orientation do not 
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per se reproduce the spatially aperiodic and complex 
structure of orientation maps in the visual cortex. 

To clarify whether the EN model is in principle cap- 
able of reproducing the biological observations, we per- 
formed an unbiased comprehensive inspection of EN 
optima for arbitrary stimulus distributions possessing 
finite fourth moments. This analysis identified two key 
parameters determining pattern selection: (i) the effec- 
tive intracortical interaction range and (ii) the fourth 
moment of the orientation stimulus distribution. We 
derived complete phase diagrams summarizing pattern 
selection in the EN model for fixed as well as variable 
retinotopy. Small interaction ranges together with low 
fourth moment values lead to either pinwheel-free 
orientation stripes, rhombic or hexagonal crystalline 
orientation map layouts as optimal states. Large interac- 
tion ranges together with orientation stimulus distribu- 
tions with high fourth moment values lead to the 
stabilization of irregular aperiodic OPM layouts. These 
solutions belong to a class of solutions previously called 
n-ECPs. This solution class encompasses a large variety 
of OPM layouts and has been identified as optimal solu- 
tions of abstract variational models of OPM develop- 
ment [35]. We showed that in the EN model due to a 
lack of a so-called permutation symmetry, among this 
family of solutions, states with low pinwheel densities 
are selected as global minima. In the extreme and pre- 
viously unexplored parameter regime of very large effec- 
tive interaction ranges and stimulus ensemble 
distributions with infinite fourth moment, permutation 
symmetry is restored and spatially aperiodic OPM lay- 
outs with higher pinwheel density are included in the 
repertoire of optimal solutions. Only in this limit, the 
repertoire of optima reproduces the recently described 
species-insensitive OPM design [38] and quantitatively 
matches experimentally observed orientation map lay- 
outs. None of these findings depend on whether the EN 
model is considered with variable or fixed retinotopy. 

Comparison to previous studies 

It is an important and long-standing question, whether 
the structure of cortical maps of variables such as stimu- 
lus orientation or receptive field position can be 
explained by a simple general principle. The concept of 
dimension reduction is a prominent candidate for such 
a principle (see e.g., [58,102] for reviews) and the quali- 
tative agreement between experimental data and pre- 
vious numerical results from dimension reduction 
models [21,42,60,62-66,68,98,102-104] can be viewed as 
evidence in favor of the dimension reduction hypothesis. 
Yet comprehensive analytical investigations of dimen- 
sion reduction problems and in particular the determi- 
nation of their optimal and nearly optimal solutions 
have been impeded by the mathematical complexity of 



these problems. For the EN algorithm applied to the 
TSP, previous analytical results established the unselec- 
tive fixed point above the first bifurcation point as well 
as the parameters at which this solution becomes 
unstable [105]. Subsequent work extended these results 
to the EN model for cortical map formation. The peri- 
odicity of solutions depending on the model parameters 
has been obtained by computing the eigenvalues of the 
Hessian matrix of the energy function [63,97,106]. Hoff- 
summer et al. [72] confirmed these results, and com- 
puted the periodicity of the emerging patterns in the 
continuous EN model formulation by linear stability 
analysis of the EN gradient descent dynamics as used in 
this study. Our results extend these findings and for the 
first time provide analytical expressions for the precise 
layout of optimal and nearly optimal dimension-redu- 
cing maps. 

In the light of the qualitative agreement between 
experimental data and numerical solutions of the EN 
model previously described, it is perhaps our most sur- 
prising result that the model's optimal dimension-redu- 
cing maps are regular periodic crystalline structures or 
pinwheel-free stripe patterns in large regions of para- 
meter space. In particular, the species-insensitive pin- 
wheel statistics observed experimentally [38] are not 
exhibited by optimal solutions of the classical EN in any 
of the previously considered parameter regimes. 

Our comparison of different numerical approaches 
indicates that the differences to previous studies are 
mainly attributable to differences in the sampling of the 
stimulus manifold in the numerical optimization proce- 
dures. In their seminal publication, Durbin and Mitchi- 
son used sets of 216 stimuli from the circular stimulus 
ensemble and applied a Gauss-Seidel procedure to 
obtain stationary configurations [21]. A similar proce- 
dure was used in [104]. Quite frequently, the number of 
stimuli used for optimization is of the same order of 
magnitude as the number of model neurons or centroids 
in feature space. This provides a relatively sparse sam- 
pling of the stimulus manifold [63-65]. Finite stimulus 
sampling effects are expected to worsen when feature 
spaces of higher dimension are considered. 

The choice of small stimulus sets in previous dimen- 
sion reduction studies was imposed mainly by the lim- 
itations of computing power. Using a parallelized 
implementation of the Cholesky method for determinis- 
tic annealing [62-65] on a multicore architecture with 2 
TB working memory, we explored the dependence of 
the obtained near optimal solutions on the sampling of 
the feature space manifold over two orders of magni- 
tude. We find that, the more stimuli are sampled, the 
closer the numerically obtained configurations resemble 
our analytical predictions. Our results on the classical 
EN model with deterministic annealing suggest that in 
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the limit of large stimulus numbers, one would perfectly 
recover our analytical results both for periodic condi- 
tions or nonperiodic boundary conditions with realistic 
system sizes. This dense stimulus sampling limit is also 
readily visible in our reproduction of the original Durbin 
and Mitchison sampling and the modification of the 
predicted map structure with stimulus number (Figure 
17). The finding that computational limitations pre- 
vented Durbin and Mitchison from obtaining the genu- 
ine predictions of their dimension reduction model 
should not be viewed as diminishing the importance of 
their contribution. The dimension reduction approach 
has played a unique and extremely productive role in 
guiding the conceptualization of cortical functional 
architecture. It has established an abstract view on corti- 
cal representations without which most of our current 
theoretical knowledge about candidate theories for corti- 
cal architectures could not have been obtained. 

Our results about optimal states of the EN for the cir- 
cular and uniform stimulus ensembles however agree 
with some prior work. In [25], the gradient descent 
dynamics of the EN model was used as a model for the 
emergence and refinement of cortical maps during 
development. Simulated visual stimulus features 
included retinotopy, orientation and eye dominance. 
The numerical procedures were similar to the one 
developed in this study. Parameters were chosen such 
that 5 4 = 4/3 and a/ A « 0.366. This study found that an 
initially large number of pinwheels decayed via pairwise 
annihilation of pinwheels with opposite topological 
charge. Our analysis predicts a stripe-like OP pattern as 
optimal solution in this regime, both in the case of a 
fixed uniform retinotopy as well as with variable retino- 
topy. In our simulations, this state is reached after an 
initial phase of symmetry breaking with the generation 
of numerous pinwheels via pairwise pinwheel annihila- 
tion. Our analytical and numerical results thus confirm, 
explain, and generalize these previous findings. 

The previous results also indicated that the inclusion 
of eye dominance in the EN model slightly slows down 
but does not stop the pinwheel annihilation process (see 
[25], Figure 3). This raises the possibility that the main 
features of our analysis of optimal solutions for the EN 
model may persist when additional feature dimensions 
are taken into account. Reichl et al. in fact observed that 
models with interacting OPM and OD maps (ODMs) 
exhibit a transition from pinwheel-free stripes to peri- 
odic pinwheel crystals similar to the transitions found in 
the EN [37] and demonstrated that this transition is a 
general feature of models with interacting OPM and 
ODMs [107]. A rigorous characterization of map struc- 
tures predicted by the simultaneous optimization of 
multiple periodic feature representations such as orien- 
tation preference and OD constitutes an important goal 



for future studies. The recent study by Reichl et al. [37] 
suggests that this issue can successfully be approached 
using concepts from the nonlinear dynamics of pattern 
formation. Finally, one recent study used the continuous 
formulation of the EN model to investigate the impact 
of postnatal cortical growth on the formation of OD 
columns in cat visual cortex [69]. Consistent with our 
results, this study also observed perfectly regular stripe- 
like patterns as stationary states in gradient descent 
simulations. The dynamics of the convergence of the 
ODMs toward the stripes was modified by including 
cortical growth into the model. However, as soon as 
growth terminated, simulated ODM layouts readily con- 
verged toward regular stripes. How cortical growth 
interacts with the formation of orientation columns is 
currently not understood and represents a further inter- 
esting topic for future studies. 

Geometric relationship between retinotopic distortions 
and OPMs 

Experimental results on the geometric relationships 
between the map of visual space and the map of orien- 
tation preference are ambiguous. Optical imaging 
experiments in cat VI suggested a systematic covaria- 
tion of inhomogeneities in the RM with singularities in 
the pattern of orientation columns in optical imaging 
experiments [96]. Regions of high gradient in the map 
of visual space preferentially appeared to overlap with 
regions of high gradient of the OPM. In ferret, however, 
it has been reported that high gradient regions of the 
map of visual space correspond to regions of low gradi- 
ent in the OPM [67]. In tree shrew VI, no local rela- 
tionships between the mapping of stimulus orientation 
and position seem to exist and the map of visual space 
appears to be ordered up to very fine scales [108]. In 
line with this, single unit recordings in cat area 17 
revealed no correlation between receptive-field position 
scatter and orientation scatter across local cell ensem- 
bles [109,110]. 

Our analysis of the EN model shows that its optimal 
states exhibit a negative correlation between the rates of 
change of orientation preference and retinotopic posi- 
tion, similar to what has been observed in the ferret 
[67]. This is expected from the principle of dimension 
reduction and in agreement with the original numerical 
results by Durbin and Mitchison [21]. However, both in 
simulations of the gradient descent dynamics and in 
deterministic annealing simulations with periodic 
boundary conditions as well as in analytically obtained 
optimal solutions, deviations from a perfectly uniform 
mapping of visual space are surprisingly weak (see Fig- 
ures 9, 10, 18, and 25 in Appendix 1). 

Deterministic annealing simulations with open non- 
periodic boundary conditions showed a substantially 
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increased magnitude of retinotopic distortions. This 
raises the possibility that different behaviors observed in 
different experiments might be at least partially related 
to the influence of boundary effects. The influence of 
boundary effects is expected to decline into the interior 
of an area, in particular for large areas as VI (see [111]). 
In the bulk of VI, we thus expect only a weak coupling 
of orientation map and retinotopic distortions according 
to the EN model. In this regime, the predictions from 
models with reduced rotational symmetry (so-called 
Shift-Twist symmetry [112]) about the coupling between 
retinotopic distortions and OPMs [33] appear to be 
more promising than the weak effects resulting from the 
coverage-continuity-compromise. Consistent with the 
measurements of Das and Gilbert [96], such models pre- 
dict small but significant positive correlations between 
the rates of change of orientation preference and retino- 
topic position [33]. Moreover, the form of the retinoto- 
pic distortions in such models is predicted to differ for 
pinwheels with positive and negative topological charge 



[92]. This interesting prediction of OPM models with 
Shift-Twist symmetry deserves to be tested by measur- 
ing the receptive field center positions around the two 
types of pinwheels with single cell resolution [12]. 

Aperiodic OPMs reflect long-range intracortical 
suppression 

Our unbiased search through the space of stimulus 
ensembles with finite fourth moment revealed the exis- 
tence of spatially aperiodic optimal solutions in the EN. 
It is important to realize that the selection of these solu- 
tions is not easily viewed as resulting from an optimal 
compromise between coverage and continuity. In fact, 
the continuity parameter in the respective parameter 
regime is so small that solutions essentially maximize 
coverage (see Figures 12, 15, and 16). Instead, this phe- 
nomenon reflects a different key factor in the stabiliza- 
tion of pinwheel-rich aperiodic layouts, namely the 
dominance of long-ranged and effectively suppressive 
interactions. This is illustrated in Figure 20 which 




Figure 20 Different patterns of evoked activity for different effective ranges of intracortical interaction in the EN model. The 

component <s z |z(x)> of the orientation map z(x) in the direction of the stimulus s z is plotted as a meshed 3D graph in a 6A x 6A patch. Color 

code and height of the projection below indicate the strength of activation. The stimulus is presented in the center of the displayed cortical 

subregion. (a) Evoked activity patterns e(x, S, z, r) for small interaction range a/A = 0.1 and weakly oriented stimulus with s z = 0.01 (left) and 

strongly oriented stimulus with s z = 4 (right). rPWCs (see upper right) are optimal in this regime, (b) Evoked activity patterns e(x, S, z, r) for large 

interaction range o~/A = 0.9 and weakly oriented stimulus with s z = 0.01 (left) and strongly oriented stimulus with s z = 4 (right). Spatially 

aperiodic 8-planforms (see upper right) are optimal in this regime. A uniform retinotopy was assumed in all cases for simplicity. 
I J 
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depicts different forms of stimulus-evoked activity pat- 
terns in the EN model. For a short-range interaction 
(Figure 20a), the activity evoked by low as well as high 
orientation energy stimuli is an almost Gaussian activity 
peak located near the stimulus position. The peak is 
shallow for low (left) and sharp for high 'orientation 
energy' (right). In the corresponding parameter regime, 
square pinwheel crystals are the optimal solution of the 
EN. For a longer range of interaction where aperiodic 
OPM layouts are the optimal states, the activity evoked 
by a single point-like stimulus is qualitatively different. 
Here, the activity pattern is extended and spans several 
hypercolumns (Figure 20b). It is weakly modulated for 
low orientation energy stimuli (left) and consists of sev- 
eral distinct peaks for high orientation energy stimuli 
(right). In this regime, neurons at a distance of several 
columns compete for activity through the normalization 
term in the EN which leads to a nonlocal and effectively 
suppressive intracortical interaction. 

It is presumably not a mere coincidence that recent stu- 
dies of abstract variational models of OPM development 
[35,38,93] mathematically identified this type of interac- 
tion as a key mechanism for stabilizing realistic OPM lay- 
outs. It has been shown that all models for OPM 
development that share the basic symmetries (i) transla- 
tional symmetry (ii) rotational symmetry (iii) shift symme- 
try and (iv) permutation symmetry and in addition are 
dominated by long-range suppressive interactions, form a 
universality class that generates maps exhibiting a univer- 
sal and realistic pinwheel statistics. In such models, sup- 
pressive long-range interactions are key to stabilizing 
irregular arrangement of pinwheels, which otherwise lar- 
gely disappear or crystalize during optimization. We have 
stressed that the EN model as considered here obeys the 
symmetries (i) to (iii). In the limit of infinite orientation 
stimulus ensemble fourth moment, permutation symmetry 
(iv) is restored. The EN can thus be tuned into the above 
universality class by sending the orientation stimulus dis- 
tribution fourth moment to infinity and choosing an expo- 
nentially small continuity parameter to realize effective 
long-range coupling. Indeed, the phase diagrams for 
abstract variational models of OPM development [35] and 
those of the EN model found here are structurally very 
similar. In both cases, a rather large orientation stripe 
phase is complemented by a cascade of phase transitions 
toward more complex, aperiodic and pinwheel-rich OPM 
layouts induced by long-range suppressive interactions. 
Using abstract variational models, it has been shown 
recently that the stabilization of regular crystalline pin- 
wheel layouts can alternatively be achieved by a strong 
coupling between the map of orientation and the map of 
eye dominance [37,107]. The structure of the phase dia- 
grams of such models however appears fundamentally dif- 
ferent from the structure of the EN phase diagrams. 



The parameter regime in which the EN model's optimal 
solutions exhibit the experimentally observed pinwheel sta- 
tistics is not at all intuitive and in our opinion questions 
the conventional interpretation of the EN model for the 
formation of cortical feature maps. Firstly, the extremely 
small continuity parameter questions the fundamental role 
of a tradeoff between coverage and continuity. We note 
that such a parameter regime is currently not accessible to 
numerical simulations. In addition, an apparently funda- 
mental property for any adequate model for OPM optimi- 
zation or development, namely a Turing-type finite 
wavelength instability of the unselective state [32], is lost in 
the limit r\ — > 0. At first sight the infinite fourth value 
required may appear reminiscent of the power-law distri- 
butions for orientation energy found in the statistics of nat- 
ural images [113,114]. However, as visualized in Figure 20b, 
the essential property of the EN model in the infinite 
fourth moment regime is the occurrence of patterns of 
activity spatially extended beyond a single hypercolumn 
representing spatially localized point-like stimuli. These 
activity patterns mediate the long-range interactions 
between distant orientation columns which in turn cause 
the stability of realistic pinwheel-rich aperiodic OPM lay- 
outs. It is obvious that spatially extended stimuli provide a 
much more plausible and realistic source of extended activ- 
ity patterns in models for visual cortical development (for 
an extended discussion see [54]). Optimization models for 
cortical maps based on the representation of more complex 
spatially extended visual stimuli, such as natural scenes, 
rather than a model based on point-like stimuli with 
extreme statistics would then be a more appropriate basis 
for understanding visual cortical functional architecture. 

Comparison to the SOFM model 

Several alternatives to the EN model have been proposed 
as optimization approaches that can account for the struc- 
ture of visual cortical maps. One prominent alternative 
dimension reduction model is the so-called SOFM, origin- 
ally introduced by Kohonen [59]. It is widely believed that 
this model, albeit lacking an exact energy functional [115], 
implements a competition between coverage and continu- 
ity very similar to the EN model [56,57,66,115]. The 
SOFM has been reported to reproduce many of the 
experimentally observed geometric properties of visual 
cortical feature maps (e.g., [56,57,61,66-68]). The numeri- 
cal procedures used in all of these studies were either the 
deterministic annealing procedure or the nonrecurring 
application of a stimulus set without systematic assess- 
ment of pattern convergence. An analysis of the nontrivial 
stationary states of a dynamical systems formulation of the 
SOFM model is currently lacking. The main difference 
between the SOFM model and the EN model is that the 
activation function by definition has the form of a stereo- 
typical Gaussian and competition is incorporated by a 



Keil and Wolf Neural Systems & Circuits 201 1, 1:17 
http://www.neuralsystemsandcircuits.eom/content/1/1/17 



Page 38 of 55 



hard winner-takes-all mechanism. As a consequence, it is 
not obvious that a long-range suppressive interaction 
regime can be realized in this model. According to our 
analysis, one would thus expect orientation stripes and 
rPWCs as nontrivial stationary states of the SOFM model. 
In a very recent study of the SOFM algorithm that used a 
numerical procedure similar to the gradient descent simu- 
lations developed in this article, both pinwheel annihila- 
tion and rhombic pinwheel crystallization have been 
observed [116]. In addition, one study that examined the 
SOFM model for orientation and retinotopy found a fast 
convergence to pinwheel-free stripe-like solutions for a 
wide parameter range [25]. In view of these results, it 
seems worthwhile to also reexamine the SOFM model 
with respect to its stationary states. 

Rugged or smooth energy landscape 

As for many optimization problems in biology, the optimi- 
zation of visual cortical functional architecture has been 
considered a problem characterized by a rugged energy 
landscape [76]. In case of the EN model the expectation of 
a rugged energy landscape at first sight seems quite plausi- 
ble. Originally, the elastic network algorithm was invented 
as a fast analogue method to approximately solve NP-hard 
problems in combinatorial optimization such as the TSP 
[77,78]. In the TSP, the stimulus positions correspond to 
the locations of cities a salesman has to visit on the short- 
est possible tour. In problems such as the TSP, the energy 
functionals to be minimized are known to possess many 
local minima and the global minimization of these func- 
tionals generally represents an extremely difficult problem 
[78]. Our analysis reveals that the trade-off between cover- 
age and continuity for the mapping of a continuous feature 
space manifold leads to a much simpler structure of the 
energy landscape. This is also indicated by the fact that 
almost all of our gradient descent dynamics simulations 
readily converged to the predicted global minimum of the 
energy functional. Figure 21 illustrates the smooth struc- 
ture of the EN energy landscape close to pattern formation 
threshold for different model parameters for a one dimen- 
sional path through the state space. In this landscape, the 
small set of stable planforms correspond to local minima 
of the EN energy functional, and unstable planforms to 
saddle points in the energy landscape. The optimal states 
correspond to global minima. Note that along the depicted 
state space path, unstable stationary solutions may appear 
as local minima if the unstable directions along which the 
energy decreases are orthogonal to the path. 

What is the origin of this qualitative difference in the 
shape of the energy landscapes? In the traveling salesman 
problem, the finite repertoire of possible tours consists of 
all permutations of the N cities that the salesman has to 
visit. By self-organized competition between the aim to 
visit all cities and the aim to minimize the path length, the 



elastic network algorithm converges to a specific ordering 
of the cities that eventually yields a very short tour. Most 
likely, the qualitative difference to the EN model for visual 
cortical map architecture originates from the transition 
from a finite number of cities to a continuum. When the 
elastic network algorithm is considered with an ensemble 
of cities (or stimuli) distributed according to a continuous 
probability density function, there is no discrete repertoire 
of tours. Both, the repertoire of tours as well as the path 
through the landscape of cities or equivalently the space of 
visual stimulus features are determined by self-organiza- 
tion. The first is generated by the symmetry breaking 
mechanism that leads to the instability of the homoge- 
neous state. The second corresponds to the selection of 
one of the many nontrivial stable steady states. 

An interesting property of the EN model dynamics 
that can be inferred from the energy landscape depicted 
in Figure 21 is the type of competition between two 
stable stationary states, where both are present in the 
system with a wall or a domain boundary between 
them. The motion of the wall or domain boundary is 
predicted to proceed in the direction that increases the 
fraction of the pattern with lower energy. An example 
of such competition can be seen in Figure 23g in 
Appendix. At t = lOOr, a small domain with an sPWC 
state is present. The area of this region is gradually 
reduced over the time course of the simulation until the 
pinwheel-free optimal state is reached. 

Are simple OPM layouts an artifact of model simplicity? 

The perfectly periodic types of stationary solutions 
(stripes, crystals) that appear to dominate the classical 
EN model for retinotopy and orientation have been 
found in other models of visual cortical layouts that are 
relatively abstract. One might therefore suspect that 
they represent a mere artifact of model simplicity. One 
conceptually appealing approach where perfectly peri- 
odic layouts have been found is wiring-length minimiza- 
tion [27]. According to this hypothesis, the structure of 
an OPM can be understood by minimizing the total 
length of dendritic and axonal processes. Maps obtained 
by stopping minimization of wire length exhibited quali- 
tatively realistic layouts (see Figure 6 in [27]). Complete 
optimization, however, leads to either stripe-like pin- 
wheel-free patterns or rPWCs, identical to the ones 
obtained in our investigation of optimal solutions of the 
EN model [27]. Similarly, stripe-like and rhombic 
optima have been found in several abstract vector-field 
approaches for OPM development [31,33,117]. 

It is ruled out by two observations, that the crystalline 
and perfectly periodic optima observed in all four opti- 
mization models, the EN model, the SOFM, the wiring- 
length minimization model, and the vector-field models 
are a mere artifact of the abstract order parameter field 
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Figure 21 Illustration of the EN energy landscape close to pattern formation threshold. The variation of the energy between states of 
ideal OS, the 2-ECP state, sPWC and possible mode configurations for 8-ECPs is shown for the case that (a) the OS state has the lowest energy 
(s 4 = 0, cr/A = 0.3), (b) the sPWC state has the lowest energy (s 4 = 0, oVA = 0.1), and (c) the most anisotropic 8-ECP has lowest energy (s 4 = 100, 
a/A = 0.8). The energy values between the state are computed from a state obtained by linear interpolation between two neighboring states 
on the x-axis. Note that not all local minima in a-c correspond to a stable fixed point of the amplitude dynamics (see text). 



description of cortical selectivity patterns that is com- 
mon to these approaches. Firstly, equally simplistic 
order parameter models for OPM development with 
long-range interactions have been shown to reproduce 
spatially irregular map layouts [35,38,93]. The occur- 
rence of periodic optimal solutions is thus not a neces- 
sity in this model class. Secondly, pinwheel 
crystallization has also been observed in detailed net- 
work models for the development of OPMs, notably in 
the first ever model for the self-organization of orienta- 
tion selectivity by von der Malsburg in 1973 [14,118]. 
Thus, on the one hand the phenomenon of pinwheel 
crystallization is thus not restricted to simple order 
parameter models and on the other hand abstract and 
mathematically relatively simple models can exhibit 
complex and biologically realistic optimal solutions. 



Map rearrangement and layout optimization 

Irrespective of the optimization principle invoked to 
describe the structure of visual cortical maps, several com- 
mon features of the resulting dynamics have been 
observed. The dynamics of optimization models usually 
starts with a phase of pattern emergence, where selectivity 
to visual features arises from an initially homogeneous 
unselective or weakly selective state. As we and others 
have shown, feature maps in these models continue to 
evolve after single cell selectivities reach mature levels. In 
fact, the phase of initial pattern emergence is typically fol- 
lowed by a prolonged phase of rearrangement of selectiv- 
ities and preferences until a stable configuration is reached 
that represents a genuine optimum. This is not an excep- 
tional type of dynamics but rather constitutes the generic 
expectation for a spatially extended system [83,84]. 
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What drives the second phase of map rearrangement? 
The initial emergence of feature selectivity is predomi- 
nantly a local process in which merely neighboring units 
interact with each other to roughly match their selectiv- 
ities. In the resulting spatial layout, selectivities are 
therefore far from being optimally arranged in space 
with respect to the global organization of selectivities on 
larger scales. Depending on the interactions incorpo- 
rated in the model, local matching processes may (i) 
effectively propagate through space optimizing the pat- 
tern over gradually increasing spatial scales or (ii) dis- 
tant sites may start to directly interact with each other 
to guide a rearrangement toward a globally optimized 
pattern after their initial emergence. 

An illustrative example is provided by the emergence 
of pinwheel-sparse orientation stripes. Qualitatively, it is 
easy to see that a pattern of orientation stripes satisfies 
the continuity constraint very well. In a stripe pattern, 
preferred orientations are constant along one direction 
in space, realizing the absolute minimum of the orienta- 
tion gradient in this direction. Reaching such a config- 
uration obviously requires to select the preferred 
orientation at widely separated sites (along the stripe 
axis) to be identical. Because initially such sites develop 
independent preferred orientations, the optimized col- 
umn layout can only emerge through a secondary rear- 
rangement process. If the dominant low energy state has 
low pinwheel density, the later phase is governed by pin- 
wheel motion and pairwise pinwheel annihilation. If this 
state is pinwheel-rich, e.g., a pinwheel crystal or an aper- 
iodic pinwheel-rich state, both pinwheel annihilation 
and pinwheel creation together with a coordinated rear- 
rangement of pinwheels are expected to occur. 

The local, essentially random processes during the initial 
emergence of a fist pattern are in principle incapable of 
directly generating an optimized layout. In fact, it has been 
established that this initial so-called symmetry breaking 
phase will in general produce a random arrangement of 
selectivities of model-insensitive statistics [25,32,36]. The 
occurrence of some form of secondary reorganization is 
thus a qualitative prediction of any optimization model, 
provided that the optimal map is not seeded by an innate 
mechanism. The results presented in this study and many 
reports demonstrate that Hebbian plasticity is capable and 
often expected to achieve such rearrangements. 

In gradient descent dynamics simulations of the EN 
model for retinotopy and stimulus orientation with con- 
ventional stimulus ensembles, pinwheel densities were 
found to be strongly time-dependent after the initial col- 
umn formation (see e.g. Figures 6, 7). In particular, the 
timescale for the establishment of full orientation selectiv- 
ity and the time needed for either annihilation of a sub- 
stantial fraction of pinwheels or their crystallization into 
periodic pinwheel crystals are in the same range of tens of 



tau. A similar time dependence of pinwheel density has 
also been observed in other models for OPM development 
with periodic optima [35,38]. Pinwheel annihilation in the 
EN can be slightly slowed down by additional features 
such as retinotopy Figures 10 and 17 or OD [25,37] but 
not by orders of magnitude. For this reason, signatures of 
the periodic optima of a developmental dynamics become 
visible at rather early simulation stages. Long-term mini- 
mization is apparently not essential to express the main 
layout features of the global minimum. 

Because the main features of the dominant optimal solu- 
tions become apparent immediately after orientation selec- 
tivity saturates it appears not easy to reproduce the 
species-independent map layout in models with periodic 
crystalline optima by pattern freezing. In our simulations 
to match even only the pinwheel density, a very precise 
timing of the freezing point would be required. There is 
currently no evidence for such a freezing mechanism in 
early development. In cats and ferrets, cortical maps for 
OD, orientation or direction arise on a timescale between 
hours and a few days (e.g., [13,119,120]). The underlying 
circuits can be rapidly modified, e.g. by deprivation experi- 
ments, even on the timescale of hours [121,122] weeks 
after full selectivity has been established. Recently, evi- 
dence for long-term visual cortical circuit reorganization 
after the emergence of feature selectivity during normal 
development has emerged in diverse systems. In mouse, 
for example, activity-dependent changes induced by nor- 
mal visual experience during the critical period, i.e., long 
after the primary emergence of orientation selectivity, 
have been shown to gradually match eye-specific inputs in 
the cortex [123]. Specifically, the data from mouse indi- 
cates that preferred orientations in the two eyes initially 
often emerge unmatched and subsequently change toward 
one binocularly matched orientation preference. Because 
preferred orientations in the two eyes initially are statisti- 
cally independent, this suggests that neurons can rotate 
their orientation preferences up to at least 45° during post- 
natal development. This is reminiscent of pairing experi- 
ments in kitten visual cortex in which Fregnac and 
coworkers induced neurons to changed their preferred 
orientation by up to 90° after pairing of a visual stimulus 
with intracortical stimulation [124,125] (see also [126]). 
Also in the cat, visual cortical orientation columns in 
visual areas VI and V2 have been found to undergo rear- 
rangement during the late phase of the critical period [41]. 
In this process, columns in mutually connected regions of 
areas VI and V2 or in retinotopically matched regions in 
left and right hemisphere areas become progressively bet- 
ter matched in size. In the same species, a systematic reor- 
ganization of OD columns during postnatal development 
has been observed [69]. Essential features of this columnar 
rearrangement are reproduced by the EN model for OD 
patterns simulated in a growing domain. 
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In view of these observations, it seems unlikely that 
aperiodic orientation maps in the visual cortex represent 
frozen transient states of a developmental dynamics whose 
attracting layouts are pinwheel crystals or pinwheel free 
states. In fact, models for the activity-dependent develop- 
ment of OPMs with aperiodic optima predict only subtle 
changes of the OPM layout during the convergence after 
the establishment of selectivity [35,38]. This might also 
explain the apparent stability of cortical maps during nor- 
mal development over short periods [119]. Further studies 
of the long-term rearrangement and stabilization of corti- 
cal functional architecture are needed to exhaustively 
characterize such processes. Given the fundamental role of 
map reorganization for any optimization theory of visual 
cortical development, chronic imaging experiments track- 
ing the spatial arrangement of feature selectivities in indi- 
vidual animals beyond the emergence of selectivity and 
through later developmental stages are expected to be 
highly informative about fundamental principles of visual 
cortical optimization. 

Conclusions 

Together with recent progress on the quantitative char- 
acterization of cortical functional architecture [38,69,93], 
this study lays the foundation for a mathematically rig- 
orous and biologically informative search for optimiza- 
tion principles that successfully explain the architecture 
of columnar contour representations in the primary 
visual cortex. A mathematically controlled and quantita- 
tively precise determination of the predictions of candi- 
date optimization principles is demanded by 
accumulating evidence indicating that geometrical fea- 
tures of visual cortical representations are biologically 
laid down with a precision in the range of a few percent 
[38,127,128]. Such data is expected to substantially 
reduce the range of candidate optimization principles 
that are consistent with biological observations. In parti- 
cular, for the principle that cortical orientation maps are 
designed to optimally compromise stimulus coverage 
and feature continuity, our analysis demonstrates that 
the classical EN model for orientation preference and 
retinotopy essentially fails at explaining the biologically 
observed architecture. Our finding that the EN model 
exhibits biologically realistic optima only in a limit in 
which point-like stimuli are represented by complex 
spatially extended activity patterns corroborates that 
large-scale interactions are essential for the stabilization 
of OPM layouts with realistic geometry [35,39,87,93]. In 
the light of these results, principles for the optimal 
representation of entire visual scenes by extended corti- 
cal activity patterns appear as promising candidates for 
future studies (see also [54]). In fact, there is recent evi- 
dence that visual cortical activity becomes progressively 
better matched to the statistics of natural stimuli but 



not to simplistic artificial stimulus ensembles [129]. We 
expect the methods developed here to facilitate a com- 
prehensive characterization of such candidate principles. 

Methods 

Expansion of EN equation 

In order to analytically calculate the approximate opti- 
mal dimension-reducing mappings in the EN model 
with fixed retinotopy, an expansion of the nonlinear EN 
OPM dynamics (Equation (3)) up to third-order around 
the unselective fixed point has to be derived. This 
expansion is briefly sketched in the following. Equation 
(3) with r(x) = 0 is of the form 

d t z(x, t) = Af x [z] + r]Az(x, t), 

where M x [z] is a nonlinear functional of z(; t), para- 
meterized by the position x. Clearly, the diffusion term 
contains no nonlinear terms in z(-> t) and therefore third 
order terms of the dynamics d t z(x, t) exclusively stem 
from third order terms of the Volterra series expansion 
of the functional N x [z\ around the fixed point z(x, t)= 
0. By the Shift symmetry (Equation (8)), only third-order 
contributions of the form N 3 [z, z, z] are allowed, i.e., 



N 3 [z,z,z] ■■ 



yd wd v- 



2 JJJ ' 8z(y)8z{w)8z{\) 
Collecting all the terms yields 

11 

N 3 [z,z,z] = J2 a i N 3[ z ]> 



z(y)z(w)z(v). 
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The coefficients aj for various orientation stimulus 
ensembles are given in 'Results' section. 

Adiabatic elimination of r(x, f) 

In order to analytically calculate the approximate optimal 
dimension-reducing mappings in the EN model with vari- 
able retinotopy, an expansion of the nonlinear EN retino- 
topy and orientation map dynamics (Equations (3, 4)) up 
to third-order around the nonselective fixed point has to 
be derived and retinotopic distortions have to be adiabati- 
cally eliminated. Both of these calculations are briefly 
sketched in the following. Equation (3) is of the form 

3 t z{x, t) = Nx[z, r] + r]Az(x, t), 

where J\f x [z, r] is a nonlinear functional of z(; t) and r 
(♦, t), parameterized by the position x. The diffusion 
term contains no nonlinear terms in z(; t) and therefore 
third order terms of the dynamics of z(x, t) exclusively 
stem from third-order terms of the Volterra series 
expansion of the functional A4[z,r] around the fixed 
point {z(x, t)= 0, r(x, t)= 0}. By the shift symmetry 
(Equation (8)), only terms in form of a cubic operator 
N3[z,z,z] and a quadratic operator Q z [r, z] are allowed 
when expanding up to third order. N3 [z, z, z] is given in 
Equation (37). Q[z, r] can be calculated via 

QZ[r,z] = // d2yd2w (^^)^ [z ' t] \^^ 
and this yields 

QZ[r ' z] = (( ' 5z ll~f 2) < x ) / ^ 2 y<r(y)^(y-x)> 

"T^ / ^ 2 MxU(y) K 2(y-x)> + J d 2 y{r(y),z(yW 2 (y-x)) 

(l 5 z| 2 ) ff 

- 36 *2 a 8 II d 2 yd 2 wz(y)(r(w),K r 3 (y-x,w-x,y-w)}, 

where (♦, •) denotes the scalar product between two 
vectors and 



K r 2 (x) 
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In complete analogy, by expanding the right hand side 
of the dynamical equation for the retinotopic distortions 
(Equation (4)) up to second-order, the vector-valued 
quadratic operator Q r [z, z] can be obtained as 



2or 2 -(|5 z | 2 ) f 2 2 
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Inserting r(x) = — L" 1 [Q r [z,z]] into Q z [r, z] and using 
the linearity of L" 1 as well as the bilinearity of both, 
Q r [z,z] and Q z [r, z], yields a sum 
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where the second equal sign is valid for (\s z \ 2 ) = 2. 
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Amplitude equations from N\[z,z,z\ 
We catalog the numerous stationary solutions of Equa- 
tion (13) following [35], by considering planforms 



,ik,x 



with an even number N of modes with amplitudes Aj 
and k ; = k c (cos(2nj/N) t sm(2nj/N )). In the vicinity of a 
finite wavelength instability, where the nonselective state 
z(x) = 0 becomes unstable with respect to a band of 
Fourier modes around a finite wave number /< c -by sym- 
metry, the dynamics of the amplitudes Aj at threshold 
has the form 

N N 

A< = Ai - A, &j I A I 2 " EWr ~ A i~ < (42) 
j=i j=i 

where /" denotes the index of the mode antiparallel 
to the mode kj = — kj-, and the coefficients 

£v = C 1 " i^iMN-^l) and 
^ = (1 — 5y — 8{-j)f(\ai — <Xj\) only depend on the 
angle \a t - 0Cj\ between mode i and The angle-depen- 
dent interaction functions g{a) and f(a) are obtained 
from Equation (13) by a multi scale expansion 
[35,83,84,87] as 



/(«) 
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(43) 



(44) 



where k 0 = /c c (l, 0) and h(a) = /c c (cos a, sin a).f{a) is 
7r-periodic, since the right hand side of Equation (44) is 
invariant with respect to the transformation h(a) — > h(a 
+ ji) = -h(a). g(a) is 27T-periodic in general. If, however, 
the nonlinearity is permutation symmetric (Equation 
(34)) it can be seen from Equation (43) that g(a) is ji- 
periodic as well. 

Stability of stationary planform solutions 

To determine the stability of fixed points of the ampli- 
tude equations (42), the eigenvalues of their stability 
matrices have to be determined. In general, for any 
fixed point A = A 0 of the dynamical system A = F(A) 
with complex-valued A and F, we have to compute 
the eigenvalues of the Hermitian 2N x 2N matrix 



/ d¥ 
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For the system of amplitude equations, we obtain 
dF- ( N \ 

j± = r8 ik - S ik I Y^gij\Aj\ 2 J - AigikAk - Ai-f ik (Ak- +A k ) 

dF- ( N \ 

^ = -AigikA k - 8i- k I J2fii A i A r I • 

Stability of a solution, or more precise intrinsic stabi- 
lity is given, if all eigenvalues of M are negative definite. 
Extrinsic stability is given, if the growth of additional 
Fourier modes is suppressed. To test whether a plan- 
form solution is extrinsically stable, we introduce a test 
mode T such that 

N 

z(x) = TA x + ^VV 

i 

with kp = (cos /3, sin P)k c . We insert this ansatz into 
Equation 15 and obtain 

N 

dtT = rT _ J^gtf - Pi) \Aj\ 2 T + 0{T 2 ) 
i 

as the dynamics of the test mode I 7 , where g{fS) is the 
angle-dependent interaction function corresponding to 
N3 [z, z, z] . For the solution T = 0 to be stable, we there- 
fore obtain the condition 

N 

i 

where we assumed k^ ^ kj, kj^ . These conditions for 

intrinsic and extrinsic stability were numerically evalu- 
ated to study the stability of n-ECPs and rPWCs. 

Coupled essentially complex planforms 

In 'Results' section, we presented a closed form expres- 
sion for the retinotopic distortions associated via Equa- 
tion (28) with stationary planform solutions of Equation 
(29). Here, we sketch how to explicitly calculate this 
representation. We start with the ansatz 



z(x) = ^ mi = k c 



(45) 



for the OPM z(x). Note that this general ansatz 
includes essentially complex planforms as well as 
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rPWCs. To simplify notation, we denote the individual 
terms in Equation (41) 

QaM] = y^ijs // dV 2 ^(y-x,y- w, w - x)z(w)z(y). 

Each of the Q f [z,z], z = 1, 2, 3 can be evaluated for 
the planform ansatz in Equation (45) and we obtain 



JV ¥ 'E^" iklX = 2 G i ' E m A Ak){kj - k fe ) sin((k ; - k k )x) 

fe,;<fe 

+ 3(A ; A fe )(kj - k fe ) cos((k ; - k fe )x)} 
JV ¥ -E^" iklX = - (20r2 ~J 5z ' 2>) E ^" (k '"" kfe)2<r2 (k; - WCA/Afc) sin((kj - k*)x) 

+ 3(A ; A fe )cos((k J -k fe )x)} 
£VV J]A fe e-^ x = - <|5zl ^ 2 £ (k ; - - k fe >- 2 " W(AjA k ) sin((k, - k fe )x) 
+ 3(A ; A k )cos((k i -k fe )x)}. 



All resulting terms are proportional to either (k, - k ; ) 
sin((k/ - k ; )x) or (k, - k ; ) cos((k, - k ; )x) i * j. These func- 
tions are longitudinal modes (see Figure 3b) which have 
been identified as eigenfunctions of the linearized 
dynamics of retinotopic distortions L r [r] with eigenvalue 

X[(|k, - k ; |) = -Ik, - kf^ + e-^-^a 2 ) 

Hence, they are eigenfunctions of L" 1 ^] with eigenva- 
lue l/A£(|kj — kj|) . Using this when inserting in Equa- 
tion (28) and setting (|s z | 2 ) = 2, we obtain expression 
(31) for the retinotopic distortions belonging to an arbi- 
trary planform. 

Phase diagrams 

To compute the regions of minimal energy shown in 
Figures 6, 10, 12, 15, and 16 as well as Figures 23, 25 in 
Appendix 1, we first computed the fixed points of Equa- 
tion (42) at each point in parameter space. For n-ECPs, 
we constructed the coupling matrix g in Equation (22) 
for all mode configurations not related by any combina- 
tion of the symmetry operations: (i) Translation: 

Aj Ajjfy, (ii) Rotation: Aj -» A J+AJ , (iii) Parity: 
Aj —> A( N _jy. By Equation (22), we then computed the 
absolute values of the corresponding amplitudes. If 
J2j=i — 0 f° r ail U a valid n-ECP fixed point of 

Equation 42 was identified. Its energy was then com- 
puted by Equation (23). For orientation stripes and 
rPWCs, the derived analytical expressions for their 
energy (Equations (18, 20)) were evaluated. To analyze 
the stability of the fixed points, the conditions for 



intrinsic and extrinsic stability (see above) were numeri- 
cally evaluated. 

Numerical procedures-gradient descent simulations 

To test our analytical calculations and explore their 
range of validity, we simulated Equations (3) and (4) on 
a 64 x 64 grid with periodic boundary conditions. Simu- 
lated systems were spatially discretized with typically 8 
grid points per expected column spacing A max of the 
orientation preference pattern (see 'Results' section) to 
achieve sufficient resolution in space. Test simulations 
with finer discretization (16 and 32 grid points per 
A max ) did not lead to different results. Progression of 
time was measured in units of the intrinsic timescale r 
(see 'Results' section) of the pattern formation process. 
The integration time step St was bounded by the rele- 
vant decay time constant of the Laplacian in Equation 
(3) around k c and by the intrinsic timescale r of the sys- 
tem, using 8t = min{l/(20?7^), r/l0}. This ensured 
good approximation to the temporally continuous 
changes of the patterns. We used an Adams-Bashforth 
scheme for the first terms on the respective r.h.s. of 
Equations (3, 4). The second terms (diffusion) were trea- 
ted by spectral integration, exhibiting unconditional 
numerical stability. The stimulus positions s r were cho- 
sen to be uniformly distributed in retinal coordinates. 
The stimulus averages in Equations (3, 4) were approxi- 
mated by choosing a random representative sample of 
N s stimuli at each integration time step, with 



, k N 0 r 2 8t 

N s = max { 10 \ 



where n corresponds to the dimensions of the feature 
space in addition to the two retinal positions (in our 
case, n = 2), T 2 = (£/A max ) 2 the squared aspect ratio of 
the simulated system in units of A 2 , s s the resolution in 
feature space, N 0 the number of stimuli we required to 
approximate the cumulative effect of the ensemble of 
stimuli within each feature space voxel s n+2 . With N 0 = 
100 and s s = 0.05, we ensured a high signal-to-noise 
ratio for all the simulations. Typical values for N s were 
between 2.5 x 10 5 and 4 x 10 6 . All simulations were 
initialized with z(x, t = 0) = 10"V 27 ^( X ) and r(x, t = 0) = 
0, where the f (x) are independent identically distributed 
random numbers uniformly in [0, 1]. Different realiza- 
tions were obtained by using different stimulus samples. 

Stimuli were drawn from different distributions, each 
with (\s z \ 2 ) = 2. We considered (i) stimuli uniformly dis- 
tributed on a ring with |s z | 2 = </2 (circular stimulus 
ensemble), (ii) stimuli uniformly distributed within a cir- 
cle {s z , \s z \ < 2} (uniform stimulus ensemble), and (iii) a 
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Gaussian stimuli ensemble with 

p Sz = 1 / (2tt) exp(— \s z \ 2 /l). In addition, we considered 
mixtures of a high-fourth moment Pearson type VII dis- 
tribution and the circular stimulus ensemble. The Pear- 
son distribution is given by 



1 



aB(m - \, \) 



|2l 



1 + 



where £(♦,♦) is the Beta function [130] 



3, and m = § + such that (|s z | 2 ) 



and 

2,(| 



s z\ ) = /or equivalently s 4 = 7- 4. 

In addition to simulations in which independent sets 
of stimuli were used for evaluating the stimulus average 
in Equations (3, 4) for every time step, we also per- 
formed simulations in which the same fixed set of N sti- 
muli was used (see 'Results' section). To determine the 
time step St in these simulations, we first calculated 



(parameters as in regular simulations) which yields the 
number of stimuli presented to the model in one intrin- 
sic time unit r in regular simulations. To subject the 
network to the same number of stimuli per intrinsic 
time scale r in fixed stimulus set simulations, the inte- 
gration time step St was in this case chosen as 



St-- 



N 

mm < — r 



1 



N 5 ' 20^ 2 ' 10 



For small N, this resulted in very small integration 
steps. For very large N, time steps were identical to the 
regular simulations. Different realizations were obtained 
by different but fixed stimulus sets. In all simulations 
with fixed stimulus sets, stimuli were drawn from the 
circular stimulus ensemble. All other numerical methods 
were chosen as in regular simulations. 

Numerical procedures - solving the EN model with 
deterministic annealing 

A large body of previous study has solved the EN 
models for various aspects of visual cortical architec- 
ture for discrete fixed sets of stimuli and using deter- 
ministic annealing. We therefore also used 
deterministic annealing with fixed discrete sets of sti- 
muli to solve the EN model for the most frequently 
used stimulus distribution. This allowed us to better 
compare our analytical and numerical results based on 
the gradient descent dynamics for a continuous stimu- 
lus with prior results. For the discrete deterministic 
annealing approach, cortical maps are described by a 



collection of M centroids {y m }^ =1 C R d that can be 
represented as a D x M matrix Y = (y 1} . . . , y M ). 
Maps forming a compromise of coverage and continu- 
ity are obtained for {x n }^ =1 C R d represented asaDx 
N matrix X = (x 1? . . . , x M ). In our case, d = 4. The 
trade-off between coverage and continuity is formu- 
lated by the energy function 



+ ^tr(Y T YS). (46) 



The matrix S determines the topology of the network 
as well as the boundary conditions and is typically 
derived from a discretized derivative based on a finite 
difference scheme or stencil (see below). For large N 
and M, the energy function in Equation (46) is equiva- 
lent to the energy functional of the continuum formula- 
tion given in Equation (1) for /3 = riN and S 
implementing the discretized Laplacian operator in two 
dimensions. 

Following [62-65,71] we minimized the EN energy 
function (Equation (46)) by an iterative deterministic 
annealing algorithm, starting with a minimization for 
large a and tracking this minimum to a small value of 
<7. As in [4], we reduced a from <7 init = 0.2 to the point 
at which the amplitude of the orientation maps saturate 
(o" « 0.03), following a = 0" init x % J where / counts the 
annealing step. This choice tracks stationary solutions of 
the EN to parameters that are very far from threshold. 
For high precision tracking of solutions, we used an 
annealing rate of % = 0.999. 

At each value of O", setting the gradient of Equation 
(46) to zero yields a nonlinear system of equations 



YA = XW 



with 



A = G + at 



S + S J 



(47) 



Here, the N x M-matrix W is given by 



1 \\xn-y m 



2 a 



w nm — 



1 \\ Xn-Ym> 

and the M x M-matrix G is 



w n 



n=l 



A is a symmetric positive-definite M x M matrix. The 
M x M matrix A is symmetric and positive-definite. 
Since both G and W depend on Y, this equation is non- 
linear in Y and has to be solved iteratively. Following 
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[62-65,71], we solved Equation (47) at each value of o 
and for each iteration directly via Cholesky-factorization. 

We implemented periodic and nonperiodic boundary 
conditions by appropriate choice of the matrix S. S 
must be positive (semi)definite for the energy to be 
bounded from below. We used the canonical 2D Lapla- 
cian stencil of order 2, to construct the M x M matrix 



/-4 + 2a 1 0 0 
1 -4 + a 1 0 
0 1 -4+al 



1 0 0 ■■• 0 
10 0 ■■• 
0 •■■ 1 ■■• 



V 1 



1 -4 + 2a 1 

1 -4+al 

1 -4 10 
1 -4 1 



v2a) 



Here, a - 0 for periodic boundary conditions and a - 
1 for nonperiodic boundary conditions. In Appendix 3, 
we also present simulation results for a fourth derivative 
stencil, in which S 2 was used for the continuity term. 
We used random stimulus positions and orientations as 
well as stimuli arranged on a grid in feature space. For 
random stimuli, positions were drawn from a uniform 
distribution in [0, 1] x [0, 1]. Orientations s z were 
drawn from the circular stimulus ensemble, with \s z \ = 
0.08 as in [65]. Stimuli from grid-like ensembles were 
distributed evenly-spaced in [0, 1] x [0, 1] and contained 
2 k evenly space orientations with \s z \ = 0.08. 

To compute the energy of pinwheel-free configuration, 
we initiated simulations with a stripe-like orientation 
preference pattern with the same typical spacing as the 
observed orientation maps and annealed from a = 0.035 
to a = 0.03. 

To enable comparison between simulations with dif- 
ferent numbers of stimuli, we scaled the continuity para- 
meter proportionally to N such that the equivalent 77 
was approximately constant. The simulated domain then 
contained a comparable number of hyper columns for 
all stimulus numbers. 

Pinwheel density from simulations 

Pinwheels locations in models were identified by the 
crossings of the zero contour lines of real and imaginary 
parts of the orientation map. Estimation of local column 
spacing A(x) was done using the wavelet analysis intro- 
duced in [127,128]. In short, an overcomplete basis of 
complex Morlet wavelets at various scales and orienta- 
tions was compared to the OPM pattern at each loca- 
tion. A(x) was estimated by the scale of the best 
matching wavelet. The mean column spacing (A(x)) x of 
a given map was then calculated from the local column 
spacing by spatial averaging. For details we refer to 
[38,127,128]. Given N pw pinwheels in a simulated 



cortical area of size L 2 , we defined the pinwheel density 
[25,38,102] 



p=N ] 



pw " 



L 2 



The pinwheel density p is a dimensionless quantity 
and depends only on the layout of orientation columns. 
The pinwheel density defined in this way is large for 
patchy and small for more band-like columnar layouts. 

Appendix 1 

The impact of nonoriented stimuli 

The main text of this article contains a complete analy- 
sis of optimal dimension-reducing mappings of the EN 
model with a circular ensemble of orientation stimuli. 
These optima are simple regular orientation stripes or 
square pinwheel crystals. The circular orientation stimu- 
lus ensemble, however, contains only stimuli with a 
fixed and finite 'orientation energy' or elongation \s z \. 
This raises the question of whether the simple nature of 
the circular stimulus ensemble might restrain the realm 
of complex dynamics in the EN model. The EN 
dynamics are expected to depend on the characteristics 
of the activity patterns evoked by the stimuli and these 
will be more diverse and complex with ensembles con- 
taining a greater diversity of stimuli. Therefore, we also 
examined the EN model in detail for a richer ensemble 
of stimuli. In this ensemble, called a uniform stimulus 
ensemble in the following, orientation stimuli are uni- 
formly distributed on the disk {s z , \s z \ < 2}, a choice 
common to many previous studies, e.g., [19,25,81]. The 
uniform ensemble in particular contains unoriented sti- 
muli with \s z \ = 0. Intuitively, the presence of these 
unoriented stimuli might be expected to fundamentally 
change the importance of pinwheels in the optimal 
OPM layout. Pinwheels' population activity is untuned 
for orientation. Pinwheel centers may therefore acquire 
a key role for the representation of unoriented stimuli. 
As such an effect should be independent of retinotopic 
distortions and to aid comparison with our previous 
results, we will again start with a fixed uniform retino- 
topy r(x) = 0. 

The linear stability properties of the unselective fixed 
point are independent of the ensemble of orientation 
stimuli «|s z | 2 ) = 2 throughout this article). The coeffi- 
cients in Equation (14), however, of course depend on 
the fourth moment of the stimulus distribution. Insert- 
ing (|5^| 4 ) = 16/3 into Equation (32), we obtain 



fl 4 = —} 

a 10 



1 

2^ 



a 2 
■ a s 



an 



1 

4tto 6 . 



a 3 = 
a<s = 
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The angle-dependent interaction functions are then 
given by 

g(a) = ± (l - 2e~^ 2 - D (l - 2.-^"-)) 

+ _L (> c V(cosa-l) _ ^ + JL,-2^ sinh 4 ( 1 / 2fe 2 ff 2 cosa ) 

/(a) = (l - ^" 2fe?ff2 (cosh(2fc 2 a 2 cos a) + 2 cosh(fc 2 a 2 cos a)) + 2<r fe ' ff2 ) 



cosh(2fc 2 a 2 cos a) - l) + ^e- 2fe ' ff2 sinh 4 (l /2fe 2 a 2 cos a) 



Both functions are depicted in Figure 22 for two dif- 
ferent values of the effective intracortical interaction 
range o7A. They qualitatively resemble the functions 
depicted in Figure 5. Figure 23 displays the phase dia- 
gram of the EN model with uniform stimulus ensemble. 
As summarized in the main part of this article, it is 
almost identical to that obtained for the circular stimu- 
lus ensemble (Figure 6). Two different optimal states are 
found, square pinwheel crystals (sPWCs) and orientation 
stripes (OSs) separated by a phase border at <7/A - 0.15. 
Both fixed points are stable for all a/A. Figure 23b-k 
demonstrates, that these analytical results are confirmed 
by direct numerical simulations of Equation (3) with r 
(x) = 0. As for the circular stimulus ensemble, we also 
tested the stability of stationary n-ECP solutions with 2 
< n < 20 by numerical evaluation of the criteria for 
intrinsic and extrinsic stability (see 'Methods' section). 
We find all n-ECPs with 2 < n < 20 intrinsically 
unstable for all interaction ranges o7A. The simple 
phase space structure furthermore apparently remains 
unchanged if we consider the model far from pattern 
formation threshold as shown in Figure 24. Simulations 
bear a close resemblance to the simulations with circu- 
lar orientation stimulus ensemble (Figure 7). Either con- 
vergence to sPWC-like patterns or patterns with large 
orientation stripe domains is observed. Again, pinwheel 
annihilation in the case of large o7A is less rapid than 
close to threshold (Figure 24a,b). The linear pinwheel- 



free zones increase their size over the time course of the 
simulations, eventually leading to a stripe pattern. For 
smaller interaction ranges o7A, the OPM layout rapidly 
converges toward a crystal-like rhombic arrangement of 
pinwheels with dislocations and pinwheel density close 
to 4. 

Figure 25 shows that taking retinotopic distortions 
into account yields an almost identical picture compared 
to the circular stimulus ensemble. For small interaction 
range o7A, the analytically predicted optimum is a quad- 
ratic pinwheel crystal with pinwheel density p = 4. For 
larger o7A, the analytically predicted optimum is an 
orientation stripe pattern with pinwheel density p = 0. 
Our results are confirmed by direct simulations of Equa- 
tions (3, 4) (Figure 25b-e). The simulation results are 
virtually indistinguishable from the circular stimulus 
ensemble. 

All together, the EN dynamic given by Equations (3, 4) 
and in particular the set of ground states of the EN 
model and their stability regions appear almost identical 
when considering either a circular or a uniform orienta- 
tion stimulus ensemble. We found two different optima 
depending on the parameter regime, orientation stripes 
for larger interaction ranges and quadratic pinwheel 
crystals for shorter interaction ranges. In addition, the 
EN dynamics appears to be unchanged by the presence 
of unoriented stimuli. 

Appendix 2 

Strength of retinotopic coupling 

In our manuscript, we have shown that retinotopic dis- 
tortions only have a weak influence on the optima of 
the EN model as well as its dynamics (see Figures 10 
and 12). Here, we quantify the influence of retinotopic 
distortions on the pattern formation process by compar- 
ing the angle-dependent interaction function for 
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Figure 23 Optimal solutions of the EN model for uniform stimulus ensemble and fixed representation of visual space, (a) At criticality, 
the phase space of this model is parameterized by either the continuity parameter r\ (blue labels) or the effective interaction range cr/A (red 
labels, see text), (b, c) OPMs (b) and their power spectra (c) in a simulation of Equation (3) with r(x) = 0, r = 0.1, a/A = 0.12 (77 = 0.57) and 
uniform stimulus ensemble, (d) Analytically predicted optimum for cr/A < 0.15 (rPWC). (e) Pinwheel density time courses for four different 
simulations (parameters as in b; gray traces, individual realizations; black trace, simulation in b; red trace, mean value), (f) Mean squared 
amplitude of the stationary pattern in simulations (parameters as in b) for different values of the control parameter r (black circles) and 
analytically predicted value (solid green line), (g, h) OPMs (g) and their power spectra (h) obtained in a simulation of Equation (3) with r(x) = 0, r 
= 0.1, a/A = 0.15 (77 = 0.41) and uniform stimulus ensemble, (i) Analytically predicted optimum for cr/A > 0.15 (orientation stripes), (j) Pinwheel 
density time courses for four different simulations (parameters as in g; gray traces, individual realizations; black trace, simulation in g; red trace, 
mean value), (k) Mean squared amplitude of the stationary pattern in simulations (parameters as in g) for different values of the control 
parameter r (black circles) and analytically predicted value (solid green line). 



retinotopic coupling g r (cc) with angle-dependent interac- 
tion function of the EN model with fixed retinotopy. 
We use the ratio 

c _ HgrQIh 
C ~ IW0II2 

as a measure to quantify the influence of retinotopic 
distortions. || ♦ || 2 denotes the 2-Norm, 

ll/Olh = f f\a)da. 
Jo 

If c is larger than one, g r (a) dominates the total interac- 
tion function g r (cc) + g(a) and retinotopic distortions may 
strongly influence the layout and stability of solutions of 
the EN model On the other hand, if c is small, the solu- 
tions and their stability properties are expected to not 
change substantially when including variable retinotopy 



into the EN model. Figure 26 displays the parameter c in 
the 5 4 -(7/A-plane for the EN model at threshold for two 
different conditions, 77 = 77^ and r\ r - 0. In the latter case, 
retinotopic distortions are expected to have the strongest 
impact. However, in both cases, c « 1, in almost all of 
the parameter space, implying little influence of retinoto- 
pic deviations. Only for small cr/A and small s 4 , c is larger 
than one. As shown in Figure 12, this leads to slight 
deformations of the stability regions for rhombs, and 
stripes in this region of parameter space but does not 
result in novel optimal solutions. 

Appendix 3 

Grid-like stimulus ensembles 

References [64,65]) performed simulations with stimuli 
distributed in regular intervals in feature space, called 
grid-like ensemble. For comparison, we also performed 
deterministic annealing simulations with grid-like 
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Figure 24 Numerical analysis of the EN dynamics with uniform orientation stimulus ensemble and fixed representation of visual space 

far from pattern formation threshold, (a) OPMs and their power spectra in a representative simulation of Equation (3) with r(x) = 0, r = 0.8, ol 

A = 0.3 (77 = 0.028) and uniform stimulus ensemble, (b) Pinwheel density time courses for four different simulations (parameters as in a; gray 

traces, individual realizations; black trace, simulation in a; red trace, mean value) (c) OPMs and their power spectra in a representative simulation 

of Equation (3) with r(x) = 0 and cr/A = 0.12 (77 = 0.57), other parameters as in a. (d) Pinwheel density time courses for four different simulations 

(parameters as in c; gray traces, individual realizations; black trace, simulation in c; red trace, mean value). 
I J 



stimulus sets of varying size with nonperiodic boundary 
conditions (see 'Methods' section). For these grid-like 
stimulus patterns, a competition between stripes and 
rhombs is observed (Figure 27a). Notably, these are the 
only two stable states identified by our analysis for the 
circular stimulus ensemble. For nonperiodic boundary 
conditions, rhombic pinwheel arrangements seem ener- 
getically favored for grid-like stimuli, almost indepen- 
dently of the size of the stimulus set. The average 
pinwheel density for N = 100 x 100 x 8 stimuli was p = 
3.4 (Figure 27b). As expected from the predominantly 
rhombic arrangement of pinwheels, NN-pinwheel dis- 
tances concentrate around half the typical column spa- 
cing and pinwheel pairs at short distances are not 
observed (Figure 27c). With these features, the maps 
obtained substantially differ from the experimentally 
observed pinwheel statistics [38]. 

The discrete EN model with fourth derivative 

In previous studies of the EN model, alternative defini- 
tions of the continuity term in the EN model have been 



explored [64]. A general continuity term for the spatially 
continuous formulation of the EN for OPM and retino- 
topy is a linear differential operator which will suppress 
the emergence of high-frequencies during the EN 
dynamics. A finite-wavelength instability is expected in 
this case, although the precise expressions for the criti- 
cal a and the typical wavelength will differ. As linear 
terms do not enter in the higher-order derivatives of the 
EN functional, changing the continuity term is not 
expected to alter the stability results obtained in this 
study. 

To numerically test this expected robustness of our 
results for the EN model with discrete fixed sets of sti- 
muli (see Figures 18 and 19), we also performed simula- 
tions using deterministic annealing with a fourth 
derivative stencil (see 'Methods' section). Figure 28 illus- 
trates that the results almost perfectly match the ones 
for the second-order derivative, considered in the main 
part of this article (Figures 18, 19 and Figure 27). 

When annealing with periodic boundary conditions, 
the solutions very much resemble our gradient descent 
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Figure 25 Phase diagram of the EN model for the joint mapping of visual space and orientation with a uniform orientation stimulus 
ensemble, (a) Regions of the 77 r -o~/A-plane in which n-ECPs or rPWCs have minimal energy, (b) Pinwheel density time courses for four different 
simulations of Equations (3, 4) with r = 0.1, o~/A = 0.13 (77 = 0.51), r\ r = 77 (grey traces, individual realizations; red trace, mean value; black trace, 
realization shown in c). (c) OPMs (upper row), their power spectra (middle row), and RMs (lower row) in a simulation of Equations (3, 4); 
parameters as in b. (d) Pinwheel density time courses for four different simulations of Equations (3, 4) with r = 0.1, cr/A = 0.3 (77 = 0.03), rj r = rj 
(grey traces, individual realizations; red trace, mean value; black trace, realization shown in e). (e) OPMs (upper row), their power spectra (middle 
row), and RMs (lower row) in a simulation of Equations (3, 4); parameters as in d. 
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Figure 26 Strength of coupling between orientation map and retinotopy in the EN model, (a) Ratio of HgrOlMlgOlh in the s 4 -o/ A-plane 
for the EN model at threshold and 77 = 77 r . (b) As a, but for 77 r = 0, i.e., strongest coupling. Note the logarithmic scaling of the colormap. 
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Figure 27 The EN model with deterministic annealing and stimuli, distributed on a grid in feature space, (a) OPMs (left) and RMs (right) 
for N = 20 x 20 x 8 (upper row), N = 50 x 50 x 8 (middle row) and N = 100 x 100 x 8 (lower row) stimuli and nonperiodic boundary 
conditions (annealing rate % = 0.999). p is the continuity parameter in the conventional definition of the EN model (see 'Methods' section, 
Equation (46)) and is scaled, such that a comparable number of columns emerges in all simulation for each N. (b) Pinwheel densities of EN 
solutions for different numbers of stimuli (annealing rate % = 0.999). Crosses mark individual simulations, red line indicates average values, (c) 
Statistics of nearest neighbor pinwheel distances for pinwheels of (upper left) arbitrary and (upper right) opposite and equal charge for 
100x100x8 stimuli and nonperiodic boundary conditions, averaged over four simulations (red curves). Black curves represent fits to the 
experimental data from [38]. Lower left: SD of pinwheel densities estimated from randomly selected regions in the OPM. Black dashed curve 
indicates SD for a two-dimensional Poisson process of equal density. 



dynamics simulations with Laplacian term. The larger 
the set of stimuli, the more stripe-like are the OPMs 
obtained (Figure 28a) and consequently pinwheel densi- 
ties decrease (Figure 28b, upper right). The exponent 
for the SD is considerably lower than for the Poisson 
process (Figure 28b, upper right). Typical NN-pinwheel 
distances concentrate around half the typical column 
spacing and in particular pinwheel pairs with short dis- 
tances lack completely (Figure 28b, lower left and right). 



For nonperiodic boundary conditions and random sti- 
muli, we found that retinotopic distortions are much 
more pronounced. They however decreased with 
increasing number of stimuli. For large stimulus num- 
bers, we observed stripe-like orientation preference 
domains which are interspersed with lattice-like pin- 
wheel arrangements (see Figure 28c), lower row, upper 
left corner of the OPM). Similarly to the periodic 
boundary conditions, short distance pinwheel pairs 
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(random stimuli) 
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Figure 28 The EN model with deterministic annealing and fourth derivative stencil, (a) OPMs (left) and RMs (right) for N = 103 (upper 
row), N = 10 4 (middle row) and N = 10 5 (lower row) stimuli, non-periodic boundary conditions, and annealing rate % = 0.999. f3 is the continuity 
parameter in the conventional definition of the EN model (see 'Methods' section, Equation (46)) and is scaled, such that a comparable number 
of columns is emerging in the simulations for each stimulus set. (b) Pinwheel densities (upper left) of EN solutions, SD of pinwheel densities 
estimated from randomly selected regions in the solutions (upper right). Crosses mark individual simulations, red line indicates average values. 
Black dashed curve indicates SD for a two-dimensional Poisson process of equal density. Statistics of nearest neighbor pinwheel distances for 
pinwheels of arbitrary (lower left) and (lower right) opposite and equal charge for 10 5 random stimuli and periodic boundary conditions, 
averaged over four simulations (red curves). Black curves represent fits to the experimental data from [38]. (c) As a but for nonperiodic boundary 
conditions, (d) As b, but for non-periodic boundary conditions, (e) OPMs (left) and RMs (right) for N = 20 x 20 x 8 (upper row), N = 50 x 50 x 8 
(middle row) and N = 100 x 100 x 8 (lower row) stimuli, nonperiodic boundary conditions, annealing rates % = 0.999. (f) As b, but for 
nonperiodic boundary conditions and grid-like stimuli. 



occur much less frequently than in the experimentally 
observed maps, indicating an increased regularity in the 
pinwheel arrangements compared to realistic OPMs 
(Figure 28d, lower left and right). This regularity also 
manifests itself in a smaller exponent of the SD com- 
pared to the Poisson process (Figure 28d). 

Simulations with grid-like stimulus as, e.g., used in 
[64,65] displayed a strong tendency toward rhombic pin- 
wheel arrangements analogous to the second derivative 
case (Figure 27e,f) 



Additional material 



Additional file 1: Rhombic pinwheel crystallization in the EN model. 

The movie shows OPMs (left) as well as their power spectrum (right). In 
the left panel, colors encode preferred orientation and brightness 
orientation selectivity. The simulation of the EN model was obtained by 
gradient descent dynamics with circular stimulus ensemble and fixed 
retinotopy. The simulation was started from the unselective fixed point z 
(x, r = 0) = 0 (parameters: r = 0.1, o/A = 0.1 (77 = 0.67)). 

Additional file 2: Pinwheel annihilation in the EN model. The movie 
shows OPMs (left) as well as their power spectrum (right). In the left 
panel, colors encode preferred orientation and brightness orientation 
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selectivity. The simulation of the EN model was obtained by gradient 
descent dynamics with circular stimulus ensemble and fixed retinotopy. 
The simulation was started from the unselective fixed point z(x, f = 0) = 
0 (parameters: r = 0.1, a/A = 0.3 (77 = 0.028)). 

Additional file 3: Convergence to fractured stripes in the EN model. 

The movie shows OPMs (left) as well as their power spectrum (right). In 
the left panel, colors encode preferred orientation and brightness 
orientation selectivity. The simulation of the EN model was obtained by 
gradient descent dynamics with fixed retinotopy. The simulation was 
started from the unselective fixed point z(x, t = 0) = 0 (parameters: r = 
0.1,o7A = 0.2 (77 = 0.2), s 4 = 6). 

Additional file 4: Hexagonal pinwheel crystallization in the EN 
model. The movie shows OPMs (left) as well as their power spectrum 
(right). In the left panel, colors encode preferred orientation and 
brightness orientation selectivity. The simulation of the EN model was 
obtained by gradient descent dynamics with fixed retinotopy. The 
simulation was started from the unselective fixed point z(x, f = 0) = 0 
(parameters: r = 0.1, o7A = 0.3 (77 = 0.028), s 4 = 8). 
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